%% Compito di "Tecniche di Ottimizzazione per il controllo" del 09/02/2009
%% Ingegneria Informatica LS
%% http://www.kida.netsons.org

%% Riporto il testo (per quello che mi ricordo) e la risoluzione del compito.
%% La soluzione e' riportata "cosi' com'e'": potrebbe contenere degli errori.
%% Intellegenti pauca!

%% ESERCIZIO 1.
%% dato il sistema con funzione di trasferimento G(s) = (s^2+0.3*s+0.7)/(s^5+13*s^4+42*s^3+137*s^2+191*s+210)
%% determinare un sistema di ordine ridotto a catena aperta utilizzando il metodo del troncamento diretto. 
%% Valutare l'errore commesso.
s = tf('s');
G = (s^2+0.3*s+0.7)/(s^5+13*s^4+42*s^3+137*s^2+191*s+210);
[sistema_bilanciato,valori_singolari] = balreal(G);
% i valori singolari pių significativi sono i primi due. L'ordine del sistema ridotto sara' quindi r=2
n = 5; %ordine del sistema di partenza
r = 2; %ordine del modello ridotto
sistema_troncato = ss(sistema_bilanciato.a(1:r,1:r),sistema_bilanciato.b(1:r),sistema_bilanciato.c(1:r),sistema_bilanciato.d);
% adesso calcolo l'errore commesso in termini di norma h infinito della funzione di trasferimento del sistema errore
G = tf(sistema_bilanciato);
Gtilde = tf(sistema_troncato);
TfErrore = G-Gtilde;
normhinf(TfErrore)

% e' possibile verificare che la somma dei valori singolari scartati moltiplicata per 2 e' minore dell'errore commesso
(valori_singolari(3)+valori_singolari(4)+valori_singolari(5))*2


%% ESERCIZIO 2.
%% dato il sistema:
A = [ 0 1; 0 0];
B = [0 1]';
C = [-1 1];
%% calcolare il regolatore ottimo per ciascuno dei seguenti pesi r ={0.01,2,100}
sistema = ss(A,B,C,0);
% verifico che il sistema sia controllabile:
rank(ctrb(sistema))
% ho utilizzato una funzione kleinman che ha signature 
% [K_ottimo,P] = kleinman(A,B,C,D,R,Q)
% dove K_ottimo č il vettore dei guadagni del regolatore ottimo
% P e' la soluzione definitiva positiva dell'equazione di Riccati associata al problema di controllo ottimo
% A,B,C,D sono le matrici di stato del sistema da retroazionare, R e Q sono le matrici peso del
% funzionale J da minimizzare
[K_ottimo,P] = kleinman(A,B,C,0,0.01,C'*C);
% verifico la stabilita' del sistema retroazionato:
eig(A-B*K_ottimo)

% adesso calcolo il regolatore per i restanti casi:
[K_ottimo,P] = kleinman(A,B,C,0,2,C'*C)
eig(A-B*K_ottimo)
[K_ottimo,P] = kleinman(A,B,C,0,100,C'*C)
eig(A-B*K_ottimo)


%% ESERCIZIO 3.
%% Scrivere un esempio di sistema passivo di ordine non inferiore a tre.

%% ESERCIZIO 4.
%% Progettare un controllore H_infinito per il sistema con matrici di stato:
A = [ 1 0.1; -10 0.4];
B1 = [1 1]'; B2=[1 1]';
C1 = [0.2 1]; C2 = [0.1 0]; C3 = [1 0]; D11 = 0;
D12 = 0;
D21 = 0;
D22 = 0.1;
D31 = 0.1;
D32 = 0;
%% Quindi confrontarne le prestazioni con quelle ottenibili con un controllore multiobiettivo
P = ltisys(A,[B1 B2],[C1; C2; C3],[D11 D12; D21 D22; D31 D32]);
[g_ottimo,h2_ottimo,K] = hinfmix(P,[1 1 1],[0 0 1 0])
compensatore = ss(K(1:2,1:2),K(1:2,3),K(3,1:2),K(3,3));
sistema_ol = ss(A,[B1 B2],[C1; C2; C3],[D11 D12; D21 D22; D31 D32]);
sistema_cl = feedback(sistema_ol,compensatore,2,3,+1);
G_cl = tf(sistema_cl);
% verifico le prestazioni ottenute:
normhinf(G_cl(1,1))
impulse(G_cl)

% adesso calcolo un controllore H_infinito/H_2:
[g_ottimo,h2_ottimo,K] = hinfmix(P,[1 1 1],[0 0 1 1])
compensatore = ss(K(1:2,1:2),K(1:2,3),K(3,1:2),K(3,3));
sistema_cl = feedback(sistema_ol,compensatore,2,3,+1);
G_cl_2 = tf(sistema_cl);
% verifico le prestazioni ottenute:
normhinf(G_cl_2(1,1))
normh2(G_cl_2(2,1))
impulse(G_cl)
impulse(G_cl_2)
impulse(G_cl)
% nel secondo caso e' possibile notare un lieve peggioramento della norma H_infinito,
% a fronte della minore escursione di energia delle variabili

%% ESERCIZIO 5.
%% Un contadino ha a disposizione due campi, uno di 300 e uno di 500 acri.
%% Puo' coltivare 6 tipi di cereali, numerati da 1 a 6.

%%  Ogni quintale di cereale necessita di una certa area (in acri) sul campo1 e sul campo2 e di una
%%  certa quantita di acqua (in metri cubi) secondo la seguente tabella (ho scritto direttamente i vincoli, individuando 12 variabili):
A = [0.02 0.04 0.01 0.06 0.02 0.02 0 0 0 0 0 0; 0 0 0 0 0 0 0.02 0.04 0.01 0.06 0.02 0.02; 200 100 50 140 115 120 200 100 50 140 115 120];
b = [ 300; 500; 400000]; %sono disponibili al massimo 400000 litri di acqua

%% Massimizzare il profitto, espresso dalla seguente funzione obiettivo:
f = [50; 60; 15; 56; 12; 23; 50; 60; 15; 56; 12; 23];

Aeq = zeros(1,12);
beq = 0;
LB = zeros(12,1);
Result=linprog(-f,A,b,Aeq,beq,LB)
% il massimo profitto si ha coltivando solo il cereale 2.