%parametri_simulazione_nuovo
rho_a = 1.2; %densità aria in kg m^-3 
R_a = 287.1; % costante dei gas per l'aria secca in J/(kg*K)
c_a = 1005; %calore specifico aria in J(kgK)^-1 
L = 2.257; % calore latente di vaporizzazione in JKg^-1
h_i = 8; % W/m^2K coefficente di convezione interna o conduttanza
h_e = 25; % W/m^2K coefficente di convezione esterna o conduttanza
s = 0.2 ; % spessore delle pareti
lambda = 0.28 ; % coefficente di conduttività (0.28=malta di calce cemento,1.40=calcestruzzo colato 0.84=muratura) W/mK
h = 2.30; l = 0.85; %altezza e larghezza porte in m
portaA = h*l; portaL =(h+h+l+l)*0.005;
g = 9.81; %accelerazione gravità
pm = 100000;%100000; %pressione media in pascal
pressioneint = pm*h *g/R_a;
K_est = 1/((1/h_i)+(1/h_e)+(s/lambda));% trasmittanza termica K per le pareti esterne in w/m^2K
K_int = 1/((2/h_i)+(s/lambda));% trasmittanza termica K per le pareti interne
K_fin =3;% 3; % trasmittanza termica per vetro camera inW/m^2K
Kter = 1.4; %1.4; %conduttività termica del terreno dipendente dalla composizione dello stesso; da varia 0.7 a 2.1
%c_mura = 1000; %calore specifico delle mura
a = 1; % coefficente di assorbimento per lunghezze d'onda corte <3 x10^-6m deimateriali da costruzione
f = 1; % fattore di trasparenza dei vetri alle radiazioni solari
finestraA1= 1*0.8; finestraL1 = (1+1+0.8+0.8)*0.001; % area e giunto finestra piccola m^2
finestraA2= 2*0.8; finestraL2 = (2+2+0.8+0.8)*0.001; % area e giunto finestra grande m^2
apri_fin= 0.0; %probabilità che le finestre siano aperte valore da 0 a 1
apri_porta= 0.0; %probabilità che le porte tra le stanze siano aperte valore da 0 a 1
%rho_mura= 1000; %densità muratura

% disturbi esterni
radiazione_solare_max =300; % Watt
T_ter=17+273.5; %273.5+10;%
T_est_max= 10+273.5;%293.5;
T_est_min= 0+273.5;%293.5
U_est_max=24*10^-3;% 24*10^-3;%%kg di h2o in kg aria secca
U_est_min= 20*10^-3; %20*10^-3;% Le due vel si sommano, se non c'è vento minimo lo
V_vento_min = 0*60/1000; %velocità vento m/min
V_vento_max = 0*60/1000; % velocità max vento in metri/minuto
D_vento = 60*4; % il vento dura 4h

% fonti interne
media_persone=2;%(abs(randn));%5;
%media_persone= floor(media_persone);
sigma_persone=2;
tv = 300; %watt
frig = 320; wash = 50;
pc = 290; pers_watt = 64; lux = 50; 
pers_vap = 60/(60*1000);%60/(60*1000); % kg h2O prodotti al minuto
doccia = 2600/(60*1000); cucina = 100/(1000*60); fiore = 7.5/(60*1000);
pianta = 11/(60*1000);


porte= 0;%istante in cui si comincia a calcolare differenza di pressione tra le stanze

%matrice di adiacenza tra le stanze
ad=[0 8*3 0
    8*3 0 8*3
    0 8*3 0];
% matrice delle prote tra le stanze
door = [0 1 0
        1 0 1
        0 1 0];
       
%parametri stanze
%pavimento, soffitto, nord, sud, est, ovest


%piano terra


%stanza 1
volume1 = 7*8*3;
cost1= rho_a*c_a*volume1;
%prima riga superfici finestre, seconda riga giunti finestre nord, sud, est, ovest
finestre1=[finestraA2 finestraA2 0 finestraA2
           finestraL2 finestraL2 0 finestraL2];
sommafinestre1=[sum(finestre1(1,:)') sum(finestre1(2,:)')];       
pareti_esterne1=[7*8 7*8 3*7-finestre1(1,1) 7*3-finestre1(1,2) 0-finestre1(1,3) 8*3-finestre1(1,4)]; %il piano terra è considerato parete esterna
B1= 8 ; b1= 7 ; % lato maggiore e minore del pavimento
costpavimento1 =((4*Kter*exp(0.5*b1/B1))/(b1*pi))*log((b1+s+b1)/(s));%0;%((4*Kter*exp(0.5*b1/B1))/(b1*pi))*atan(b1/(b1+s));
somma1=sum(pareti_esterne1')-pareti_esterne1(1,1);
pareti_interne1=ad(1,:);
[riga1,colonna1,area1]=find(pareti_interne1);


%stanza 2
volume2 = 9*8*3;
cost2= rho_a*c_a*volume2;
%superfici finestre nord, sud, est, ovest
finestre2=[2*finestraA2 2*finestraA2 0 0
           2*finestraL2 2*finestraL2 0 0];
sommafinestre2=[sum(finestre2(1,:)') sum(finestre2(2,:)')];         
pareti_esterne2=[9*8 9*8 3*9-finestre2(1,1) 9*3-finestre2(1,2) 0-finestre2(1,3) 0-finestre2(1,4)]; %il piano terra è considerato parete esterna
B2= 9 ; b2= 8 ; % lato maggiore e minore del pavimento
costpavimento2 =((4*Kter*exp(0.5*b2/B2))/(b2*pi))*log((b2+b2+s)/(s));
somma2=sum(pareti_esterne2')-pareti_esterne2(1,1);
pareti_interne2=ad(2,:);
[riga2,colonna2,area2]=find(pareti_interne2);


%stanza 3
volume3 = 4*8*3;%4*8*3
cost3= rho_a*c_a*volume3;
%superfici finestre nord, sud, est, ovest
finestre3=[0 finestraA2 finestraA2 0
           0 finestraL2 finestraL2 0];
sommafinestre3=[sum(finestre3(1,:)') sum(finestre3(2,:)')];         
pareti_esterne3=[4*8 4*8 4*3-finestre3(1,1) 4*3-finestre3(1,2) 8*3-finestre3(1,3) 0-finestre3(1,4)]; %il piano terra è considerato parete esterna
B3= 8 ; b3= 4 ; % lato maggiore e minore del pavimento
costpavimento3 =((4*Kter*exp(0.5*b3/B3))/(b3*pi))*log((b3+b3+s)/(s));
somma3=sum(pareti_esterne3')-pareti_esterne3(1,1);
pareti_interne3=ad(3,:);
[riga3,colonna3,area3]=find(pareti_interne3);




% calcolo dei punti di equilibrio e matrici per la linearizzazione
% fonti interne per il calcolo dell'ingresso di equilibrio
fonti_interne1 = [0,0];%[5*pers_watt+pc+4*lux+tv 5*pers_vap+5*fiore+pianta];%[0,0];%
fonti_interne2 = [0,0];%[5*pers_watt+5*pc+6*lux 5*pers_vap+5*fiore+pianta];
fonti_interne3 = [0,0];%[5*pers_watt+1*pc+4*lux+tv 5*pers_vap+1*fiore+5*pianta];



Veq = [(radiazione_solare_max/2) (T_est_max+T_est_min)/2 (U_est_max+U_est_min)/2 (V_vento_max+V_vento_min) T_ter]';
Xeq =[273.5+22 273.5+22 273.5+22 0.023 0.023 0.023]';% [22+273.5 22+273.5 22+273.5 0.024 0.024 0.024]'; Xeq = complex(Xeq);
Ueq = [000 0 000 0 000 0 ]'; Ueq = complex (Ueq);

Ueq(2,1) = -(cost1/c_a)*((sommafinestre1(1,2)/volume1)*(Veq(3,1)-Xeq(4,1))*sign(-(1/Veq(2,1))+(1/Xeq(1,1)))*(pressioneint*abs((-1/Veq(2,1))+(1/Xeq(1,1))))^(3/2)+(c_a/cost1)*fonti_interne1(1,2)+(1/volume1)*(Veq(3,1)-Xeq(4,1))*(finestre1(2,1)*(0.5*rho_a*Veq(4,1)^2)^(3/2)+finestre1(2,2)*(0.5*rho_a*(Veq(4,1)/2)^2)^(3/2))+(1/volume1)*door(1,2)*portaL*(Xeq(5,1)-Xeq(4,1))*sign(-(1/Xeq(2,1))+(1/Xeq(1,1)))*(pressioneint*(abs(-(1/Xeq(2,1))+(1/Xeq(1,1)))))^(3/2));
Ueq(1,1) = L*Ueq(2,1)-cost1*((sommafinestre1(1,2)/volume1)*(Veq(2,1)-Xeq(1,1))*sign(-(1/Veq(2,1))+(1/Xeq(1,1)))*(pressioneint*abs(-(1/Veq(2,1))+(1/Xeq(1,1))))^(3/2)+((K_int*ad(1,2))/cost1)*(Xeq(2,1)-Xeq(1,1))+costpavimento1*(pareti_esterne1(1,1)/cost1)*(Veq(5,1)-Xeq(1,1))+(K_fin*sommafinestre1(1,1)/cost1)*(Veq(2,1)-Xeq(1,1))+(K_est*somma1/cost1)*(Veq(2,1)-Xeq(1,1))+(K_est/cost1)*(pareti_esterne1(1,2)+pareti_esterne1(1,4))*(Veq(1,1)*a/h_e)+(finestre1(1,2))*(f*Veq(1,1))/cost1+(finestre1(2,1)*(0.5*rho_a*Veq(4,1)^2)^(3/2)+finestre1(2,2)*(0.5*rho_a*(Veq(4,1)/2)^2)^(3/2))*(Veq(2,1)-Xeq(1,1))/volume1+(door(1,2)*portaL*(pressioneint^(3/2))/volume1)*(Xeq(2,1)-Xeq(1,1))*sign(-(1/Xeq(2,1))+(1/Xeq(1,1)))*(abs(-(1/Xeq(2,1))+(1/Xeq(1,1))))^(3/2)+(1/cost1)*fonti_interne1(1,1));
Ueq(4,1) = -(cost2/c_a)*(((sommafinestre2(1,2)/volume2)*(Veq(3,1)-Xeq(5,1))*sign(-(1/Veq(2,1))+(1/Xeq(2,1)))*(pressioneint*abs(-(1/Veq(2,1))+(1/Xeq(2,1))))^(3/2)+(c_a/cost2)*fonti_interne2(1,2)+(1/volume2)*(Veq(3,1)-Xeq(5,1))*(finestre2(2,1)*(0.5*rho_a*Veq(4,1)^2)^(3/2)+finestre2(2,2)*(0.5*rho_a*(Veq(4,1)/2)^2)^(3/2))+(1/volume2)*door(2,1)*portaL*(Xeq(4,1)-Xeq(5,1))*sign(-(1/Xeq(1,1))+(1/Xeq(2,1)))*(pressioneint*(abs(-(1/Xeq(1,1))+(1/Xeq(2,1)))))^(3/2)+(1/volume2)*door(2,3)*portaL*(Xeq(6,1)-Xeq(5,1))*sign(-(1/Xeq(3,1))+(1/Xeq(2,1)))*(pressioneint*(abs(-(1/Xeq(3,1))+(1/Xeq(2,1)))))^(3/2)));
Ueq(3,1) = L*Ueq(4,1)-cost2*((sommafinestre2(1,2)/volume2)*(Veq(2,1)-Xeq(2,1))*sign(-(1/Veq(2,1))+(1/Xeq(2,1)))*(pressioneint*abs(-(1/Veq(2,1))+(1/Xeq(2,1))))^(3/2)+(K_int*ad(2,1)/cost2)*(Xeq(1,1)-Xeq(2,1))+(K_int*ad(2,3)/cost2)*(Xeq(3,1)-Xeq(2,1))+costpavimento2*(pareti_esterne2(1,1)/cost2)*(Veq(5,1)-Xeq(2,1))+(K_fin*sommafinestre2(1,1)/cost2)*(Veq(2,1)-Xeq(2,1))+(K_est*somma2/cost2)*(Veq(2,1)-Xeq(2,1))+(K_est/cost2)*(pareti_esterne2(1,2)+pareti_esterne2(1,4))*(Veq(1,1)*a/h_e)+(finestre2(1,2))*f*Veq(1,1)/cost2+(finestre2(2,1)*(0.5*rho_a*Veq(4,1)^2)^(3/2)+finestre2(2,2)*(0.5*rho_a*(Veq(4,1)/2)^2)^(3/2))*(Veq(2,1)-Xeq(2,1))/volume2+(door(2,1)*portaL*pressioneint^(3/2)/volume2)*(Xeq(1,1)-Xeq(2,1))*sign(-(1/Xeq(1,1))+(1/Xeq(2,1)))*(abs(-(1/Xeq(1,1))+(1/Xeq(2,1))))^(3/2)+(door(2,3)*portaL*pressioneint^(3/2)/volume2)*(Xeq(3,1)-Xeq(2,1))*sign(-(1/Xeq(3,1))+(1/Xeq(2,1)))*(abs(-(1/Xeq(3,1))+(1/Xeq(2,1))))^(3/2)+(1/cost2)*fonti_interne2(1,1));
Ueq(6,1) = -(cost3/c_a)*((sommafinestre3(1,2)/volume3)*(Veq(3,1)-Xeq(6,1))*sign(-(1/Veq(2,1))+(1/Xeq(3,1)))*(pressioneint*abs(-(1/Veq(2,1))+(1/Xeq(3,1))))^(3/2)+(c_a/cost3)*fonti_interne3(1,2)+(1/volume3)*(Veq(3,1)-Xeq(6,1))*(finestre3(2,1)*(0.5*rho_a*Veq(4,1)^2)^(3/2)+finestre3(2,2)*(0.5*rho_a*(Veq(4,1)/2)^2)^(3/2))+(1/volume3)*door(3,2)*portaL*(Xeq(5,1)-Xeq(6,1))*sign(-(1/Xeq(2,1))+(1/Xeq(3,1)))*(pressioneint*(abs(-(1/Xeq(2,1))+(1/Xeq(3,1)))))^(3/2));
Ueq(5,1) = L*Ueq(6,1)-cost3*((sommafinestre3(1,2)/volume3)*(Veq(2,1)-Xeq(3,1))*sign(-(1/Veq(2,1))+(1/Xeq(3,1)))*(pressioneint*abs(-(1/Veq(2,1))+(1/Xeq(3,1))))^(3/2)+((K_int*ad(3,2))/cost3)*(Xeq(2,1)-Xeq(3,1))+costpavimento3*(pareti_esterne3(1,1)/cost3)*(Veq(5,1)-Xeq(3,1))+(K_fin*sommafinestre3(1,1)/cost3)*(Veq(2,1)-Xeq(3,1))+(K_est*somma3/cost3)*(Veq(2,1)-Xeq(3,1))+(K_est/cost3)*(pareti_esterne3(1,2)+pareti_esterne3(1,4))*(Veq(1,1)*a/h_e)+(finestre3(1,2))*(f*Veq(1,1))/cost3+(finestre3(2,1)*(0.5*rho_a*Veq(4,1)^2)^(3/2)+finestre3(2,2)*(0.5*rho_a*(Veq(4,1)/2)^2)^(3/2))*(Veq(2,1)-Xeq(3,1))/volume3+(door(3,2)*portaL*(pressioneint^(3/2))/volume3)*(Xeq(2,1)-Xeq(3,1))*sign(-(1/Xeq(2,1))+(1/Xeq(3,1)))*(abs(-(1/Xeq(2,1))+(1/Xeq(3,1))))^(3/2)+(1/cost3)*fonti_interne3(1,1));


% Matrice A

format long e;
A=[(sommafinestre1(1,2)*(-1)*(pressioneint)^(3/2)/volume1)*(-sign(-1/Xeq(1,1))*abs(-1/Xeq(1,1))^(3/2)+(-Xeq(1,1))*(3/2)*(Xeq(1,1))^(-2)*abs(-1/Xeq(1,1))^(1/2))-((K_int*sum(area1))/cost1)-costpavimento1*pareti_esterne1(1,1)/cost1-K_fin*sommafinestre1(1,1)/cost1-K_est*somma1/cost1+door(1,2)*(1/volume1)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(2,1)-1/Xeq(1,1))*(abs(1/Xeq(2,1)-1/Xeq(1,1))^(3/2))+(Xeq(2,1)-Xeq(1,1))*((3/2)*(Xeq(1,1)^(-2))*((abs(1/Xeq(2,1)-1/Xeq(1,1)))^(1/2))))+door(1,3)*(1/volume1)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(3,1)-1/Xeq(1,1))*(abs(1/Xeq(3,1)-1/Xeq(1,1)))^(3/2)+(Xeq(3,1)-Xeq(1,1))*(3/2)*((abs(1/Xeq(3,1)-1/Xeq(1,1)))^(1/2))*Xeq(1,1)^(-2)) K_int*pareti_interne1(1,2)/cost1+door(1,2)*(1/volume1)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(2,1)-1/Xeq(1,1))*abs(1/Xeq(2,1)-1/Xeq(1,1))^(3/2)-(Xeq(2,1)-Xeq(1,1))*((3/2)*(abs(1/Xeq(2,1)-1/Xeq(1,1))^(1/2))*Xeq(2,1)^(-2))) K_int*pareti_interne1(1,3)/cost1+door(1,3)*(1/volume1)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(3,1)-1/Xeq(1,1))*abs(1/Xeq(3,1)-1/Xeq(1,1))^(3/2)-(Xeq(3,1)-Xeq(1,1))*((3/2)*(abs(1/Xeq(3,1)-1/Xeq(1,1))^(1/2))*Xeq(3,1)^(-2))) 0 0 0
  K_int*pareti_interne2(1,1)/cost2+door(2,1)*(1/volume2)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(1,1)-1/Xeq(2,1))*abs(1/Xeq(1,1)-1/Xeq(2,1))^(3/2)-(Xeq(1,1)-Xeq(2,1))*((3/2)*(abs(1/Xeq(1,1)-1/Xeq(2,1)))^(1/2)*Xeq(1,1)^(-2))) (sommafinestre2(1,2)*(-1)*(pressioneint)^(3/2)/volume2)*(-sign(-1/Xeq(2,1))*abs(-1/Xeq(2,1))^(3/2)+(-Xeq(2,1))*(3/2)*(Xeq(2,1))^(-2)*abs(-1/Xeq(2,1))^(1/2))-((K_int*sum(area2))/cost2)-costpavimento2*pareti_esterne2(1,1)/cost2-K_fin*sommafinestre2(1,1)/cost2-K_est*somma2/cost2+door(2,1)*(1/volume2)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(1,1)-1/Xeq(2,1))*(abs(1/Xeq(1,1)-1/Xeq(2,1)))^(3/2)+(Xeq(1,1)-Xeq(2,1))*(3/2)*((abs(1/Xeq(1,1)-1/Xeq(2,1)))^(1/2))*(Xeq(2,1)^(-2)))+door(2,3)*(1/volume2)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(3,1)-1/Xeq(2,1))*(abs(1/Xeq(3,1)-1/Xeq(2,1)))^(3/2)+(Xeq(3,1)-Xeq(2,1))*((3/2)*(abs(1/Xeq(3,1)-1/Xeq(2,1))^(1/2))*Xeq(2,1)^(-2))) K_int*pareti_interne2(1,3)/cost2+door(2,3)*(1/volume2)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(3,1)-1/Xeq(2,1))*abs(1/Xeq(3,1)-1/Xeq(2,1))^(3/2)-(Xeq(3,1)-Xeq(2,1))*((3/2)*(abs(1/Xeq(3,1)-1/Xeq(2,1))^(1/2))*(Xeq(3,1)^(-2)))) 0 0 0
  K_int*pareti_interne3(1,1)/cost3+door(3,1)*(1/volume3)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(1,1)-1/Xeq(3,1))*abs(1/Xeq(1,1)-1/Xeq(3,1))^(3/2)-(Xeq(1,1)-Xeq(3,1))*((3/2)*(abs(1/Xeq(1,1)-1/Xeq(3,1))^(1/2))*Xeq(1,1)^(-2))) K_int*pareti_interne3(1,2)/cost3+door(3,2)*(1/volume3)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(2,1)-1/Xeq(3,1))*abs(1/Xeq(2,1)-1/Xeq(3,1))^(3/2)-(Xeq(2,1)-Xeq(3,1))*((3/2)*(abs(1/Xeq(2,1)-1/Xeq(3,1))^(1/2))*Xeq(2,1)^(-2))) (sommafinestre3(1,2)*(-1)*(pressioneint)^(3/2)/volume3)*(-sign(-1/Xeq(3,1))*abs(-1/Xeq(3,1))^(3/2)+(-Xeq(3,1))*(3/2)*(Xeq(3,1))^(-2)*abs(-1/Xeq(3,1))^(1/2))-((K_int*sum(area3))/cost3)-costpavimento3*pareti_esterne3(1,1)/cost3-K_fin*sommafinestre3(1,1)/cost3-K_est*somma3/cost3+door(3,1)*(1/volume3)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(1,1)-1/Xeq(3,1))*abs(1/Xeq(1,1)-1/Xeq(3,1))^(3/2)+(Xeq(1,1)-Xeq(3,1))*((3/2)*(abs(1/Xeq(1,1)-1/Xeq(3,1))^(1/2))*Xeq(3,1)^(-2)))+door(3,2)*(1/volume3)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(2,1)-1/Xeq(3,1))*abs(1/Xeq(2,1)-1/Xeq(3,1))^(3/2)+(Xeq(2,1)-Xeq(3,1))*((3/2)*(abs(1/Xeq(2,1)-1/Xeq(3,1))^(1/2))*Xeq(3,1)^(-2))) 0 0 0 
  (sommafinestre1(1,2)*(-1)*(pressioneint)^(3/2)/volume1)*(-Xeq(4,1))*((3/2)*(Xeq(1,1))^(-2)*abs(-1/Xeq(1,1))^(1/2))+(1/volume1)*door(1,2)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(5,1)-Xeq(4,1)))*(sign(1/Xeq(2,1)-1/Xeq(1,1))*3/2*(abs(1/Xeq(2,1)-1/Xeq(1,1))^(1/2))*(Xeq(1,1)^(-2)))+(1/volume1)*door(1,3)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(6,1)-Xeq(4,1)))*(sign(1/Xeq(3,1)-1/Xeq(1,1))*3/2*(abs(1/Xeq(3,1)-1/Xeq(1,1)))^(1/2)*(Xeq(1,1)^(-2))) (1/volume1)*door(1,2)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(5,1)-Xeq(4,1)))*(sign(1/Xeq(2,1)-1/Xeq(1,1))*3/2*(abs(1/Xeq(2,1)-1/Xeq(1,1)))^(1/2)*(-Xeq(2,1)^(-2))) (1/volume1)*door(1,3)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(6,1)-Xeq(4,1)))*(sign(1/Xeq(3,1)-1/Xeq(1,1))*3/2*(abs(1/Xeq(3,1)-1/Xeq(1,1)))^(1/2)*(-(Xeq(3,1))^(-2))) (sommafinestre1(1,2)*(-1)*(pressioneint)^(3/2)/volume1)*(-sign(-(1/Xeq(1,1)))*abs(-(1/Xeq(1,1)))^(3/2))-(1/volume1)*door(1,2)*portaL*sign(1/Xeq(2,1)-1/Xeq(1,1))*((-1)*(pressioneint)*abs(1/Xeq(2,1)-1/Xeq(1,1)))^(3/2)-(1/volume1)*door(1,3)*portaL*sign(1/Xeq(3,1)-1/Xeq(1,1))*((-1)*(pressioneint)*abs(-1/Xeq(3,1)+1/Xeq(1,1)))^(3/2) (1/volume1)*door(1,2)*portaL*sign(1/Xeq(2,1)-1/Xeq(1,1))*((-1)*(pressioneint)*abs(1/Xeq(2,1)-1/Xeq(1,1)))^(3/2) (1/volume1)*door(1,3)*portaL*sign(1/Xeq(3,1)-1/Xeq(1,1))*((-1)*(pressioneint)*abs(1/Xeq(3,1)-1/Xeq(1,1)))^(3/2)
  (1/volume2)*door(2,1)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(4,1)-Xeq(5,1)))*(sign(1/Xeq(1,1)-1/Xeq(2,1))*3/2*(abs(1/Xeq(1,1)-1/Xeq(2,1)))^(1/2)*(-Xeq(1,1)^(-2))) (sommafinestre2(1,2)*(-1)*(pressioneint)^(3/2)/volume2)*(-Xeq(5,1))*((3/2)*(Xeq(2,1))^(-2)*abs(-1/Xeq(2,1))^(1/2))+(1/volume2)*door(2,1)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(4,1)-Xeq(5,1)))*(sign(1/Xeq(1,1)-1/Xeq(2,1))*3/2*(abs(1/Xeq(1,1)-1/Xeq(2,1))^(1/2))*(Xeq(2,1)^(-2)))+(1/volume2)*door(2,3)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(6,1)-Xeq(5,1)))*(sign(1/Xeq(3,1)-1/Xeq(2,1))*3/2*(abs(1/Xeq(3,1)-1/Xeq(2,1)))^(1/2)*(Xeq(2,1)^(-2))) (1/volume2)*door(2,3)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(6,1)-Xeq(5,1)))*(sign(1/Xeq(3,1)-1/Xeq(2,1))*3/2*(abs(1/Xeq(3,1)-1/Xeq(2,1)))^(1/2)*(-(Xeq(3,1))^(-2))) (1/volume2)*door(2,1)*portaL*sign(1/Xeq(1,1)-1/Xeq(2,1))*((-1)*(pressioneint)*abs(1/Xeq(1,1)-1/Xeq(2,1)))^(3/2) (sommafinestre2(1,2)*(-1)*(pressioneint)^(3/2)/volume2)*(-sign(-(1/Xeq(2,1)))*abs(-(1/Xeq(2,1)))^(3/2))-(1/volume2)*door(2,1)*portaL*sign(1/Xeq(1,1)-1/Xeq(2,1))*((-1)*(pressioneint)*abs(1/Xeq(1,1)-1/Xeq(2,1)))^(3/2)-(1/volume2)*door(2,3)*portaL*sign(1/Xeq(3,1)-1/Xeq(2,1))*((-1)*(pressioneint)*abs(1/Xeq(3,1)-1/Xeq(2,1)))^(3/2) (1/volume2)*door(2,3)*portaL*sign(1/Xeq(3,1)-1/Xeq(2,1))*((-1)*(pressioneint)*abs(1/Xeq(3,1)-1/Xeq(2,1)))^(3/2)
  (1/volume3)*door(3,1)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(4,1)-Xeq(6,1)))*(3/2*sign(1/Xeq(1,1)-1/Xeq(3,1))*abs(1/Xeq(1,1)-1/Xeq(3,1))^(1/2)*(-(1/Xeq(1,1))^(1/2))) (1/volume3)*door(3,2)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(5,1)-Xeq(6,1)))*(3/2*sign(1/Xeq(2,1)-1/Xeq(3,1))*abs(1/Xeq(2,1)-1/Xeq(3,1))^(1/2)*(-(1/Xeq(2,1))^(1/2))) (sommafinestre3(1,2)*(-1)*(pressioneint)^(3/2)/volume3)*(-Xeq(6,1))*((3/2)*(Xeq(3,1))^(-2)*abs(-1/Xeq(3,1))^(1/2))+(1/volume3)*door(3,1)*portaL*(((-1)*(pressioneint)^(3/2))*(Xeq(4,1)-Xeq(6,1)))*(3/2*sign(1/Xeq(1,1)-1/Xeq(3,1))*abs(1/Xeq(1,1)-1/Xeq(3,1))^(1/2)*((Xeq(3,1)^(-2))))+(1/volume3)*door(3,2)*portaL*(((-1)*(pressioneint)^(3/2))*(Xeq(5,1)-Xeq(6,1)))*(3/2*sign(1/Xeq(2,1)-1/Xeq(3,1))*abs(1/Xeq(2,1)-1/Xeq(3,1))^(1/2)*((Xeq(3,1)^(-2)))) (1/volume3)*door(3,1)*portaL*sign(1/Xeq(1,1)-1/Xeq(3,1))*((-1)*(pressioneint)*abs(1/Xeq(1,1)-1/Xeq(3,1)))^(3/2) (1/volume3)*door(3,2)*portaL*sign(1/Xeq(2,1)-1/Xeq(3,1))*((-1)*(pressioneint)*abs(1/Xeq(2,1)-1/Xeq(3,1)))^(3/2) (sommafinestre3(1,2)*(-1)*(pressioneint)^(3/2)/volume3)*(-sign(-(1/Xeq(3,1)))*abs(-(1/Xeq(3,1)))^(3/2))-(1/volume3)*door(3,1)*portaL*sign(1/Xeq(1,1)-1/Xeq(3,1))*((-1)*(pressioneint)*abs(1/Xeq(1,1)-1/Xeq(3,1)))^(3/2)-(1/volume3)*door(3,2)*portaL*sign(1/Xeq(2,1)-1/Xeq(3,1))*((-1)*(pressioneint)*abs(1/Xeq(2,1)-1/Xeq(3,1)))^(3/2)];

 
A1=[(sommafinestre1(1,2)*(-1)*(pressioneint)^(3/2)/volume1)*(-sign((1/Veq(2,1))-1/Xeq(1,1))*abs((1/Veq(2,1))-1/Xeq(1,1))^(3/2)+(-Xeq(1,1))*(3/2)*(Xeq(1,1))^(-2)*abs((1/Veq(2,1))-1/Xeq(1,1))^(1/2))-((K_int*sum(area1))/cost1)-costpavimento1*pareti_esterne1(1,1)/cost1-K_fin*sommafinestre1(1,1)/cost1-K_est*somma1/cost1+door(1,2)*(1/volume1)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(2,1)-1/Xeq(1,1))*(abs(1/Xeq(2,1)-1/Xeq(1,1))^(3/2))+(Xeq(2,1)-Xeq(1,1))*((3/2)*(Xeq(1,1)^(-2))*((abs(1/Xeq(2,1)-1/Xeq(1,1)))^(1/2))))+door(1,3)*(1/volume1)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(3,1)-1/Xeq(1,1))*(abs(1/Xeq(3,1)-1/Xeq(1,1)))^(3/2)+(Xeq(3,1)-Xeq(1,1))*(3/2)*((abs(1/Xeq(3,1)-1/Xeq(1,1)))^(1/2))*Xeq(1,1)^(-2)) K_int*pareti_interne1(1,2)/cost1+door(1,2)*(1/volume1)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(2,1)-1/Xeq(1,1))*abs(1/Xeq(2,1)-1/Xeq(1,1))^(3/2)-(Xeq(2,1)-Xeq(1,1))*((3/2)*(abs(1/Xeq(2,1)-1/Xeq(1,1))^(1/2))*Xeq(2,1)^(-2))) K_int*pareti_interne1(1,3)/cost1+door(1,3)*(1/volume1)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(3,1)-1/Xeq(1,1))*abs(1/Xeq(3,1)-1/Xeq(1,1))^(3/2)-(Xeq(3,1)-Xeq(1,1))*((3/2)*(abs(1/Xeq(3,1)-1/Xeq(1,1))^(1/2))*Xeq(3,1)^(-2))) 0 0 0
  K_int*pareti_interne2(1,1)/cost2+door(2,1)*(1/volume2)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(1,1)-1/Xeq(2,1))*abs(1/Xeq(1,1)-1/Xeq(2,1))^(3/2)-(Xeq(1,1)-Xeq(2,1))*((3/2)*(abs(1/Xeq(1,1)-1/Xeq(2,1)))^(1/2)*Xeq(1,1)^(-2))) (sommafinestre2(1,2)*(-1)*(pressioneint)^(3/2)/volume2)*(-sign((1/Veq(2,1))-1/Xeq(2,1))*abs((1/Veq(2,1))-1/Xeq(2,1))^(3/2)+(-Xeq(2,1))*(3/2)*(Xeq(2,1))^(-2)*abs((1/Veq(2,1))-1/Xeq(2,1))^(1/2))-((K_int*sum(area2))/cost2)-costpavimento2*pareti_esterne2(1,1)/cost2-K_fin*sommafinestre2(1,1)/cost2-K_est*somma2/cost2+door(2,1)*(1/volume2)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(1,1)-1/Xeq(2,1))*(abs(1/Xeq(1,1)-1/Xeq(2,1)))^(3/2)+(Xeq(1,1)-Xeq(2,1))*(3/2)*((abs(1/Xeq(1,1)-1/Xeq(2,1)))^(1/2))*(Xeq(2,1)^(-2)))+door(2,3)*(1/volume2)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(3,1)-1/Xeq(2,1))*(abs(1/Xeq(3,1)-1/Xeq(2,1)))^(3/2)+(Xeq(3,1)-Xeq(2,1))*((3/2)*(abs(1/Xeq(3,1)-1/Xeq(2,1))^(1/2))*Xeq(2,1)^(-2))) K_int*pareti_interne2(1,3)/cost2+door(2,3)*(1/volume2)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(3,1)-1/Xeq(2,1))*abs(1/Xeq(3,1)-1/Xeq(2,1))^(3/2)-(Xeq(3,1)-Xeq(2,1))*((3/2)*(abs(1/Xeq(3,1)-1/Xeq(2,1))^(1/2))*(Xeq(3,1)^(-2)))) 0 0 0
  K_int*pareti_interne3(1,1)/cost3+door(3,1)*(1/volume3)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(1,1)-1/Xeq(3,1))*abs(1/Xeq(1,1)-1/Xeq(3,1))^(3/2)-(Xeq(1,1)-Xeq(3,1))*((3/2)*(abs(1/Xeq(1,1)-1/Xeq(3,1))^(1/2))*Xeq(1,1)^(-2))) K_int*pareti_interne3(1,2)/cost3+door(3,2)*(1/volume3)*portaL*((-1)*(pressioneint)^(3/2))*(sign(1/Xeq(2,1)-1/Xeq(3,1))*abs(1/Xeq(2,1)-1/Xeq(3,1))^(3/2)-(Xeq(2,1)-Xeq(3,1))*((3/2)*(abs(1/Xeq(2,1)-1/Xeq(3,1))^(1/2))*Xeq(2,1)^(-2))) (sommafinestre3(1,2)*(-1)*(pressioneint)^(3/2)/volume3)*(-sign((1/Veq(2,1))-1/Xeq(3,1))*abs((1/Veq(2,1))-1/Xeq(3,1))^(3/2)+(-Xeq(3,1))*(3/2)*(Xeq(3,1))^(-2)*abs((1/Veq(2,1))-1/Xeq(3,1))^(1/2))-((K_int*sum(area3))/cost3)-costpavimento3*pareti_esterne3(1,1)/cost3-K_fin*sommafinestre3(1,1)/cost3-K_est*somma3/cost3+door(3,1)*(1/volume3)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(1,1)-1/Xeq(3,1))*abs(1/Xeq(1,1)-1/Xeq(3,1))^(3/2)+(Xeq(1,1)-Xeq(3,1))*((3/2)*(abs(1/Xeq(1,1)-1/Xeq(3,1))^(1/2))*Xeq(3,1)^(-2)))+door(3,2)*(1/volume3)*portaL*((-1)*(pressioneint)^(3/2))*(-sign(1/Xeq(2,1)-1/Xeq(3,1))*abs(1/Xeq(2,1)-1/Xeq(3,1))^(3/2)+(Xeq(2,1)-Xeq(3,1))*((3/2)*(abs(1/Xeq(2,1)-1/Xeq(3,1))^(1/2))*Xeq(3,1)^(-2))) 0 0 0 
  (sommafinestre1(1,2)*(-1)*(pressioneint)^(3/2)/volume1)*(Veq(3,1)-Xeq(4,1))*((3/2)*(Xeq(1,1))^(-2)*abs((1/Veq(2,1))-1/Xeq(1,1))^(1/2))+(1/volume1)*door(1,2)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(5,1)-Xeq(4,1)))*(sign(1/Xeq(2,1)-1/Xeq(1,1))*3/2*(abs(1/Xeq(2,1)-1/Xeq(1,1))^(1/2))*(Xeq(1,1)^(-2)))+(1/volume1)*door(1,3)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(6,1)-Xeq(4,1)))*(sign(1/Xeq(3,1)-1/Xeq(1,1))*3/2*(abs(1/Xeq(3,1)-1/Xeq(1,1)))^(1/2)*(Xeq(1,1)^(-2))) (1/volume1)*door(1,2)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(5,1)-Xeq(4,1)))*(sign(1/Xeq(2,1)-1/Xeq(1,1))*3/2*(abs(1/Xeq(2,1)-1/Xeq(1,1)))^(1/2)*(-Xeq(2,1)^(-2))) (1/volume1)*door(1,3)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(6,1)-Xeq(4,1)))*(sign(1/Xeq(3,1)-1/Xeq(1,1))*3/2*(abs(1/Xeq(3,1)-1/Xeq(1,1)))^(1/2)*(-(Xeq(3,1))^(-2))) (sommafinestre1(1,2)*(-1)*(pressioneint)^(3/2)/volume1)*(-sign((1/Veq(2,1))-(1/Xeq(1,1)))*abs((1/Veq(2,1))-(1/Xeq(1,1)))^(3/2))-(1/volume1)*door(1,2)*portaL*sign(1/Xeq(2,1)-1/Xeq(1,1))*((-1)*(pressioneint)*abs(1/Xeq(2,1)-1/Xeq(1,1)))^(3/2)-(1/volume1)*door(1,3)*portaL*sign(1/Xeq(3,1)-1/Xeq(1,1))*((-1)*(pressioneint)*abs(-1/Xeq(3,1)+1/Xeq(1,1)))^(3/2) (1/volume1)*door(1,2)*portaL*sign(1/Xeq(2,1)-1/Xeq(1,1))*((-1)*(pressioneint)*abs(1/Xeq(2,1)-1/Xeq(1,1)))^(3/2) (1/volume1)*door(1,3)*portaL*sign(1/Xeq(3,1)-1/Xeq(1,1))*((-1)*(pressioneint)*abs(1/Xeq(3,1)-1/Xeq(1,1)))^(3/2)
  (1/volume2)*door(2,1)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(4,1)-Xeq(5,1)))*(sign(1/Xeq(1,1)-1/Xeq(2,1))*3/2*(abs(1/Xeq(1,1)-1/Xeq(2,1)))^(1/2)*(-Xeq(1,1)^(-2))) (sommafinestre2(1,2)*(-1)*(pressioneint)^(3/2)/volume2)*(Veq(3,1)-Xeq(5,1))*((3/2)*(Xeq(2,1))^(-2)*abs((1/Veq(2,1))-1/Xeq(2,1))^(1/2))+(1/volume2)*door(2,1)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(4,1)-Xeq(5,1)))*(sign(1/Xeq(1,1)-1/Xeq(2,1))*3/2*(abs(1/Xeq(1,1)-1/Xeq(2,1))^(1/2))*(Xeq(2,1)^(-2)))+(1/volume2)*door(2,3)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(6,1)-Xeq(5,1)))*(sign(1/Xeq(3,1)-1/Xeq(2,1))*3/2*(abs(1/Xeq(3,1)-1/Xeq(2,1)))^(1/2)*(Xeq(2,1)^(-2))) (1/volume2)*door(2,3)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(6,1)-Xeq(5,1)))*(sign(1/Xeq(3,1)-1/Xeq(2,1))*3/2*(abs(1/Xeq(3,1)-1/Xeq(2,1)))^(1/2)*(-(Xeq(3,1))^(-2))) (1/volume2)*door(2,1)*portaL*sign(1/Xeq(1,1)-1/Xeq(2,1))*((-1)*(pressioneint)*abs(1/Xeq(1,1)-1/Xeq(2,1)))^(3/2) (sommafinestre2(1,2)*(-1)*(pressioneint)^(3/2)/volume2)*(-sign((1/Veq(2,1))-(1/Xeq(2,1)))*abs((1/Veq(2,1))-(1/Xeq(2,1)))^(3/2))-(1/volume2)*door(2,1)*portaL*sign(1/Xeq(1,1)-1/Xeq(2,1))*((-1)*(pressioneint)*abs(1/Xeq(1,1)-1/Xeq(2,1)))^(3/2)-(1/volume2)*door(2,3)*portaL*sign(1/Xeq(3,1)-1/Xeq(2,1))*((-1)*(pressioneint)*abs(1/Xeq(3,1)-1/Xeq(2,1)))^(3/2) (1/volume2)*door(2,3)*portaL*sign(1/Xeq(3,1)-1/Xeq(2,1))*((-1)*(pressioneint)*abs(1/Xeq(3,1)-1/Xeq(2,1)))^(3/2)
  (1/volume3)*door(3,1)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(4,1)-Xeq(6,1)))*(3/2*sign(1/Xeq(1,1)-1/Xeq(3,1))*abs(1/Xeq(1,1)-1/Xeq(3,1))^(1/2)*(-(1/Xeq(1,1))^(1/2))) (1/volume3)*door(3,2)*portaL*(((-1)*(pressioneint)^3/2)*(Xeq(5,1)-Xeq(6,1)))*(3/2*sign(1/Xeq(2,1)-1/Xeq(3,1))*abs(1/Xeq(2,1)-1/Xeq(3,1))^(1/2)*(-(1/Xeq(2,1))^(1/2))) (sommafinestre3(1,2)*(-1)*(pressioneint)^(3/2)/volume3)*(Veq(3,1)-Xeq(6,1))*((3/2)*(Xeq(3,1))^(-2)*abs((1/Veq(2,1))-1/Xeq(3,1))^(1/2))+(1/volume3)*door(3,1)*portaL*(((-1)*(pressioneint)^(3/2))*(Xeq(4,1)-Xeq(6,1)))*(3/2*sign(1/Xeq(1,1)-1/Xeq(3,1))*abs(1/Xeq(1,1)-1/Xeq(3,1))^(1/2)*((Xeq(3,1)^(-2))))+(1/volume3)*door(3,2)*portaL*(((-1)*(pressioneint)^(3/2))*(Xeq(5,1)-Xeq(6,1)))*(3/2*sign(1/Xeq(2,1)-1/Xeq(3,1))*abs(1/Xeq(2,1)-1/Xeq(3,1))^(1/2)*((Xeq(3,1)^(-2)))) (1/volume3)*door(3,1)*portaL*sign(1/Xeq(1,1)-1/Xeq(3,1))*((-1)*(pressioneint)*abs(1/Xeq(1,1)-1/Xeq(3,1)))^(3/2) (1/volume3)*door(3,2)*portaL*sign(1/Xeq(2,1)-1/Xeq(3,1))*((-1)*(pressioneint)*abs(1/Xeq(2,1)-1/Xeq(3,1)))^(3/2) (sommafinestre3(1,2)*(-1)*(pressioneint)^(3/2)/volume3)*(-sign((1/Veq(2,1))-(1/Xeq(3,1)))*abs((1/Veq(2,1))-(1/Xeq(3,1)))^(3/2))-(1/volume3)*door(3,1)*portaL*sign(1/Xeq(1,1)-1/Xeq(3,1))*((-1)*(pressioneint)*abs(1/Xeq(1,1)-1/Xeq(3,1)))^(3/2)-(1/volume3)*door(3,2)*portaL*sign(1/Xeq(2,1)-1/Xeq(3,1))*((-1)*(pressioneint)*abs(1/Xeq(2,1)-1/Xeq(3,1)))^(3/2)];
  
  
B=[1/cost1 -L/cost1 0 0 0 0
   0 0 1/cost2 -L/cost2 0 0
   0 0 0 0 1/cost3 -L/cost3
   0 c_a/cost1 0 0 0 0
   0 0 0 c_a/cost2 0 0 
   0 0 0 0 0 c_a/cost3 ];

C= eye(6);
  
raggiungibilita= ctrb(A,B);
raggiungibilita=rank(raggiungibilita);
osservabilita = obsv (A,C);
osservabilita = rank(osservabilita);
eigenvalue = eig (A);

alfa=[1 1 1 1 1 1]'; % coefficienti di importanza relativadelle di delta delle variabili controllate: se si vuole che le variazioni medie quadratiche di deltax2 sia 10 volte più piccola in m.q. di deltax1 allora alfa2 = 10^2
beta = [1 1 1 1 1 1]'; % coefficienti di importanza reladiva di delta delle variabili di controllo
%deltax= [abs(Xeq(1,1)-288.8) abs(Xeq(2,1)-286.8) abs(Xeq(3,1)-287.25) abs(Xeq(4,1)-0.0227) abs(Xeq(5,1)-0.0226) abs(Xeq(6,1)-0.0228)]'; % variazione massima  stimata in modulo per le variabili controllate rispetto al punto di equilibrio;
deltax= [3 3 3 0.002 0.002 0.002]'; % variazione massima  stimata in modulo per le variabili controllate rispetto al punto di equilibrio;
deltau = [2000 0.001 2000 0.001 2000 0.001  ]';% variazione massima  stimata in modulo per il controllo rispetto a quello di equilibrio;
%deltau = [abs(Ueq(1,1)+674.8642) abs(Ueq(2,1)+0.00081) abs(Ueq(3,1)+1746.3566) abs(Ueq(4,1)+0.00077) abs(Ueq(5,1)+822.6827) abs(Ueq(6,1)+0.000104) ]';% variazione massima  stimata in modulo per il controllo rispetto a quello di equilibrio;


Q = 1*eye (6);
%Q = 100*Q;
Q = diag([alfa(1,1)/deltax(1,1)^2 alfa(2,1)/deltax(2,1)^2 alfa(3,1)/deltax(3,1)^2 alfa(4,1)/deltax(4,1)^2 alfa(5,1)/deltax(5,1)^2 alfa(6,1)/deltax(6,1)^2]);
R =1*eye (6);
%R =R*100;
R = diag([beta(1,1)/deltau(1,1)^2 beta(2,1)/deltau(2,1)^2 beta(3,1)/deltau(3,1)^2 beta(4,1)/deltau(4,1)^2 beta(5,1)/deltau(5,1)^2 beta(6,1)/deltau(6,1)^2 ]);
[G,K,e] = lqr(A,B,Q,R); % K soluzione dell'equazione di riccati, G è la matrice dei guadagni di kalman
eigenchiusa= eig(A-B*G);
