%clear all %close all %clc %% Terms for the matrices in the state space model % Simplification terms: m_C_tilde = 2*c_C*c_H*m_sf_bar/(n_H*(2*c_C*m_sf_bar/n_C+2))-2*c_C^2*c_H*m_sf_bar^2/(n_C*n_H*(2*c_C*m_sf_bar/n_C+2)^2); m_g_tilde = (n_CO2-2*n_O)*m_C_tilde/n_C+(n_H2O-n_O)*c_H/(2*n_H)+eta_h2o+c_sf_O; m_C_bar = c_C*c_H*m_sf_bar^2/(n_H*(2*c_C*m_sf_bar/n_C+2))+chi*gamma_1*m_a_bar+gamma_2*M_C_char_bar; m_g_bar = (n_CO2-2*n_O)*m_C_bar/n_C+((n_H2O-n_O)*c_H/(2*n_H)+eta_h2o+c_sf_O)*m_sf_bar+(c_N+c_a_O)*m_a_bar; % Water equation: tw1 = M_w*c_w; tw2 = m_w*c_iw; tw3 = 4*epsilon_A_bcw*sigma*T_bc_bar^3+alpha_A_bcw; tw4 = alpha_A_gcw; tw5 = alpha_A_egw; tw6 = 4*epsilon_A_bcw*sigma*T_w_bar^3+alpha_A_bcw+alpha_A_gcw+alpha_A_egw+m_w*c_w; % Burning chamber equation: tbc1 = M_bc*c_bc+M_cp*c_cp; tbc3 = 4*epsilon_A_bcw*sigma*T_bc_bar^3+alpha_A_bcw+m_g_bar*c_bc; tbc6 = 4*epsilon_A_bcw*sigma*T_w_bar^3+alpha_A_bcw; tbc7 = m_C_tilde*((h_c_1+3*h_c_2)/4)+m_C_tilde*n_CO*h_c_5/(4*n_C)+c_H*h_c_6-m_g_tilde*c_bc*T_bc_bar; tbc8 = c_a*T_a+((h_c_1+3*h_c_2)/4+n_CO*h_c_5/(4*n_C))*chi*gamma_1-((n_CO2-2*n_O)*chi*gamma_1/n_C+c_N+c_a_O)*c_bc*T_bc_bar; tbc9 = ((h_c_1+3*h_c_2)/4+n_CO*h_c_5/(4*n_C)-(n_CO2-2*n_O)*c_bc*T_bc_bar/n_C)*gamma_2; % Gas chamber equation: tgc1 = M_gc*c_gc; tgc3 = m_g_bar*c_bc; tgc4 = m_g_bar*c_gc+alpha_A_gcw; tgc6 = alpha_A_gcw; tgc7 = (c_bc*T_bc_bar-c_gc*T_gc_bar)*m_g_tilde; tgc8 = ((n_CO2-2*n_O)*chi*gamma_1/n_C+c_N+c_a_O)*(c_bc*T_bc_bar-c_gc*T_gc_bar); tgc9 = ((n_CO2-2*n_O)*gamma_2/n_C)*(c_bc*T_bc_bar-c_gc*T_gc_bar); % Exhaust gas equation: teg1 = M_eg*c_eg; teg4 = m_g_bar*c_gc; teg5 = m_g_bar*c_eg+alpha_A_egw; teg6 = alpha_A_egw; teg7 = (c_gc*T_gc_bar-c_eg*T_eg_bar)*m_g_tilde; teg8 = ((n_CO2-2*n_O)*chi*gamma_1/n_C+c_N+c_a_O)*(c_gc*T_gc_bar-c_eg*T_eg_bar); teg9 = ((n_CO2-2*n_O)*gamma_2/n_C)*(c_gc*T_gc_bar-c_eg*T_eg_bar); % Solid fuel equation: mchar7 = (-m_C_tilde+c_C); mchar8 = chi*gamma_1; mchar9 = gamma_2; % O2 equation: o27 = c_sf_O-2*n_O*m_C_tilde/n_C-c_H*n_O/(2*n_H)-O2_bar*m_g_tilde; o28 = 2*c_a_O-2*n_O*chi*gamma_1/n_C-((n_CO2-2*n_O)*chi*gamma_1/n_C+c_N+c_a_O)*O2_bar; o29 = (n_O+(n_CO2-2*n_O)*O2_bar)*gamma_2/n_C; o210 = m_g_bar; %% State space model M = [tw1 0 0 0 0 0 0 tbc1 0 0 0 0 0 0 tgc1 0 0 0 0 0 0 teg1 0 0 0 0 0 0 1 0 0 0 0 0 0 M_total]; A_tilde = [-tw6 tw3 tw4 tw5 0 0 tbc6 -tbc3 0 0 tbc9 0 tgc6 tgc3 -tgc4 0 tgc9 0 teg6 0 teg4 -teg5 teg9 0 0 0 0 0 -mchar9 0 0 0 0 0 -o29 -o210]; A = inv(M)*A_tilde; B_tilde = [0 0 tw2 tbc7 tbc8 0 tgc7 tgc8 0 teg7 teg8 0 mchar7 -mchar8 0 o27 o28 0]; B = inv(M)*B_tilde; B_new =[B(:,1) B(:,2)]; Bi = [B_new zeros(2)]; E = [B(:,3)]; Ei = [E 0 0]; C = [1 0 0 0 0 0 0 0 0 0 0 158]; Ai = [A [0 0 ; 0 0 ; 0 0 ; 0 0 ; 0 0 ; 0 0] C zeros(2)]; Ci = [C zeros(2) [0 0 0 0 0 0 ; 0 0 0 0 0 0] eye(2)]; D = [0 0 0 0 0 0]; Di = [0 0 0 0 0 0 0 0]; G_tilde = [1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 1]; G = inv(M)*G_tilde; sysSS = ss(A,B,C,D); [Ad Bd Cd Dd] = c2dm(A,B,C,D,1,'zoh'); Ad=[Ad zeros(6) zeros(6) [0; 0; 0; 0; 0; 0] [0; 0; 0; 0; 0; 0] [0; 0; 0; 0; 0; 0]; 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0]; Bd = [Bd;0 0 0 ; 0 0 0;0 0 0 ; 0 0 0;0 0 0 ; 0 0 0;0 0 0 ; 0 0 0;0 0 0 ; 0 0 0;0 0 0 ; 0 0 0;0 0 0 ; 0 0 0;0 0 0]; Cd=[[Cd(1:5) 0;0 0 0 0 0 0] [zeros(2)] [zeros(2)] [zeros(2)] [zeros(2)] [zeros(2)] [zeros(2)] [zeros(2)] [0;158]]; [Adi Bdi Cdi Ddi] = c2dm(Ai,Bi,Ci,Di,1,'zoh'); Adi = [Ad [0 0 ; 0 0 ; 0 0 ; 0 0 ; 0 0 ; 0 0;0 0 ; 0 0 ; 0 0 ; 0 0 ; 0 0 ; 0 0;0 0 ; 0 0 ; 0 0 ; 0 0 ; 0 0 ; 0 0;0 0 ; 0 0 ; 0 0] Cd eye(2)]; Cdi = [Cd zeros(2) [0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] eye(2)]; Md=[M zeros(6) zeros(6) [0; 0; 0; 0; 0; 0] [0; 0; 0; 0; 0; 0] [0; 0; 0; 0; 0; 0];0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1]; Gd = inv(Md)*eye(21); Bd_new = [Bd(:,1) Bd(:,2)]; Bdi = [Bd_new zeros(2)]; Ed = [Bd(:,3)]; Edi = [Ed 0 0]; Ddi = [0 0 0 0 0 0 0 0]; sysSSi = ss(Ai,Bi,Ci,Di); save state_space_mdl_wheat.mat %% Optimal control Qn = [1 0 0 1]; Rn = [1 0 ; 0 1]; Ld = dlqe(Ad,Gd,Cd,Cd'*Qn*Cd,Rn); Ti = 600; Qri = [1*10^1/(100^2) 0 0 0 0 1*10^7/(100^2) 0 0 0 0 (1*10^1/(100^2))/(Ti^2) 0 0 0 0 (1*10^7/(100^2))/(Ti^2)]; R = [1*10^1/(100^2) 0 ; 0 1*10^7/(100^2)]; Cdiw = Cdi; Cdiw(2,6)=158; Cdiw(2,21)=0; Fde = dlqr(Adi,Bdi,Cdiw'*Qri*Cdiw,R); Fd1 = [Fde(:,1) Fde(:,2) Fde(:,3) Fde(:,4) Fde(:,5) Fde(:,6) Fde(:,7) Fde(:,8) Fde(:,9) Fde(:,10) Fde(:,11) Fde(:,12) Fde(:,13) Fde(:,14) Fde(:,15) Fde(:,16) Fde(:,17) Fde(:,18) Fde(:,19) Fde(:,20) Fde(:,21)]; Fdi = [Fde(:,22) Fde(:,23)]; %% Operation values u_bar = [m_sf_bar m_a_bar]'; y_bar = [T_w_bar O2_bar*100]; d_bar = [T_iw_bar]; x_bar = [T_w_bar T_bc_bar T_gc_bar T_eg_bar M_C_char_bar O2_bar]; r_bar = [T_w_bar O2_bar]'; r = [73.5 10.1]'; disp('State space model for corn loaded')