%% reset all and load constants clear all clc %% load wood pills constants load('.\..\Matlab\constants_wood_pills.mat') %% 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 for wood pills M_wood = [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_wood = [-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_wood = inv(M_wood)*A_tilde_wood; B_tilde_wood = [0 0 tw2 tbc7 tbc8 0 tgc7 tgc8 0 teg7 teg8 0 mchar7 -mchar8 0 o27 o28 0]; B_wood = inv(M_wood)*B_tilde_wood; B_new_wood =[B_wood(:,1) B_wood(:,2)]; E_wood = [B_wood(:,3)]; % C_wood = [1 0 0 0 0 0 % 0 0 0 1 0 0 % 0 0 0 0 0 158]; C_wood = [1 0 0 0 0 0 0 0 0 0 0 158]; D_wood = [0 0 0 0 0 0]; G_tilde_wood = [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_wood)*G_tilde_wood; %sysSS = ss(A,B,C,D); %eig(A); % [Ad Bd Cd Dd] = c2dm(A,B,C,D,100,'zoh'); % eig(Ad); %% bar values wood u_bar_wood = [m_sf_bar m_a_bar]'; y_bar_wood = [T_w_bar O2_bar*100];%[T_w_bar T_eg_bar O2_bar*100];% d_bar_wood = [T_iw_bar]; x_bar_wood = [T_w_bar T_bc_bar T_gc_bar T_eg_bar M_C_char_bar O2_bar]; r_bar_wood = [T_w_bar O2_bar]'; save mdl_wood.mat %% Erase workspace, to make wheat model. %clear all %% load wheat pills constants load('.\..\Matlab\constants_wheat.mat') %% 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 for wheat pills M_wheat = [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_wheat = [-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_wheat = inv(M_wheat)*A_tilde_wheat; B_tilde_wheat = [0 0 tw2 tbc7 tbc8 0 tgc7 tgc8 0 teg7 teg8 0 mchar7 -mchar8 0 o27 o28 0]; B_wheat = inv(M_wheat)*B_tilde_wheat; B_new_wheat =[B_wheat(:,1) B_wheat(:,2)]; E_wheat = [B_wheat(:,3)]; % C_wheat = [1 0 0 0 0 0 % 0 0 0 1 0 0 % 0 0 0 0 0 158]; % C_wheat = [1 0 0 0 0 0 0 0 0 0 0 158]; D_wheat = [0 0 0 0 0 0]; G_tilde_wheat = [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_wheat)*G_tilde_wheat; %% bar values wheat u_bar_wheat = [m_sf_bar m_a_bar]'; y_bar_wheat = [T_w_bar O2_bar*100]; %T_eg_bar d_bar_wheat = [T_iw_bar]; x_bar_wheat = [T_w_bar T_bc_bar T_gc_bar T_eg_bar M_C_char_bar O2_bar]; r_bar_wheat = [T_w_bar O2_bar]'; %% Data % % Step wood pills: load('.\..\stoker\xPC_Data_20071217_174202.mat') %Sim. time: 8500 m_sf = [xpcdata6.time(5001:90000)-500 xpcdata6.Stok(5001:90000)]; m_a = [xpcdata6.time(5001:90000)-500 xpcdata6.Blow(5001:90000)]; T_w = [xpcdata6.time(5001:90000)-500 xpcdata6.Temp(5001:90000)]; T_iw = [xpcdata6.time(5001:90000)-500 xpcdata6.Retur(5001:90000)]; T_eg = [xpcdata6.time(5001:90000)-500 xpcdata6.Tr0xF8g(5001:90000)]; O2 = [xpcdata6.time(5001:90000)-500 xpcdata6.Lambda(5001:90000)]; fan = [xpcdata6.time(5001:90000)-500 xpcdata6.Fan(5001:90000)]; %%%set beta motor value to wood pills. beta_motor_corn = 117.07*10^-6; beta_motor_wood = 95.4789499218371*10^-6; %%%set input data and output u_corn = [xpcdata6.time(5001:90000)-500 xpcdata6.Stok(5001:90000)*beta_motor_corn xpcdata6.Blow(5001:90000)*beta_blower]; u_wood = [xpcdata6.time(5001:90000)-500 xpcdata6.Stok(5001:90000)*beta_motor_wood xpcdata6.Blow(5001:90000)*beta_blower]; d = [xpcdata6.time(5001:90000)-500 xpcdata6.Retur(5001:90000)]; y_meas = [xpcdata6.time(5001:90000)-500 xpcdata6.Temp(5001:90000) xpcdata6.Lambda(5001:90000)]; %xpcdata6.Tr0xF8g(5001:90000) %% Steady state wood pills % load('.\..\stoker\xPC_Data_20071217_121801.mat') %Sim. time: 4450 % m_sf = [xpcdata6.time(12001:56845)-xpcdata6.time(12001) xpcdata6.Stok(12001:56845)]; % m_a = [xpcdata6.time(12001:56845)-xpcdata6.time(12001) xpcdata6.Blow(12001:56845)]; % T_w = [xpcdata6.time(12001:56845)-xpcdata6.time(12001) xpcdata6.Temp(12001:56845)]; % T_iw = [xpcdata6.time(12001:56845)-xpcdata6.time(12001) xpcdata6.Retur(12001:56845)]; % T_eg = [xpcdata6.time(12001:56845)-xpcdata6.time(12001) xpcdata6.Tr0xF8g(12001:56845)]; % O2 = [xpcdata6.time(12001:56845)-xpcdata6.time(12001) xpcdata6.Lambda(12001:56845)]; % fan = [xpcdata6.time(12001:56845)-xpcdata6.time(12001) xpcdata6.Fan(12001:56845)]; % % set beta motor value to wood pills. % beta_motor = 95.4789499218371*10^-6; % % set input data and output % u = [xpcdata6.time(12001:56845)-xpcdata6.time(12001) xpcdata6.Stok(12001:56845)*beta_motor xpcdata6.Blow(12001:56845)*beta_blower]; % d = [xpcdata6.time(12001:56845)-xpcdata6.time(12001) xpcdata6.Retur(12001:56845)]; % y_meas = [xpcdata6.time(12001:56845)-xpcdata6.time(12001) xpcdata6.Temp(12001:56845) xpcdata6.Lambda(12001:56845)]; xpcdata6.Tr0xF8g(12001:56845) %% Load data (step corn) % load('.\..\stoker\xPC_Data_20080502_163308.mat') % sim time: 18946 % m_sf = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Stok(102500:end)]; % m_a = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Blow(102500:end)]; % T_w = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Temp(102500:end)]; % T_iw = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Retur(102500:end)]; % T_eg = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Tr0xF8g(102500:end)]; % O2 = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Lambda(102500:end)]; % fan = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Fan(102500:end)]; % % %set beta motor value to corn. % beta_motor_corn = 117.07*10^-6; % beta_motor_wood = 95.4789499218371*10^-6; % %set input data and output % % u_corn = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Stok(102500:end)*beta_motor_corn xpcdata6.Blow(102500:end)*beta_blower]; % u_wood = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Stok(102500:end)*beta_motor_wood xpcdata6.Blow(102500:end)*beta_blower]; % d = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Retur(102500:end)]; % y_meas = [xpcdata6.time(102500:end)-xpcdata6.time(102500) xpcdata6.Temp(102500:end) xpcdata6.Lambda(102500:end)]; % xpcdata6.Tr0xF8g(102500:end) %% load data (mixed fuels) % Mixed fuels. % sim time: 4465 % load('.\..\stoker\xPC_Data_20080513_154324.mat') %Sim. time: 4450 % m_sf = [xpcdata6.time xpcdata6.Stok]; % m_a = [xpcdata6.time xpcdata6.Blow]; % T_w = [xpcdata6.time xpcdata6.Temp]; % T_iw = [xpcdata6.time xpcdata6.Retur]; % T_eg = [xpcdata6.time xpcdata6.Tr0xF8g]; % O2 = [xpcdata6.time xpcdata6.Lambda]; % fan = [xpcdata6.time xpcdata6.Fan]; % % %set beta motor value to corn. % beta_motor_wheat = 117.07*10^-6; % beta_motor_corn = 117.07*10^-6; % beta_motor_wood = 95.4789499218371*10^-6; % %set input data and output % % u_corn = [xpcdata6.time xpcdata6.Stok*beta_motor_corn xpcdata6.Blow*beta_blower]; % u_wood = [xpcdata6.time xpcdata6.Stok*beta_motor_wood xpcdata6.Blow*beta_blower]; % d = [xpcdata6.time xpcdata6.Retur]; % y_meas = [xpcdata6.time xpcdata6.Temp xpcdata6.Lambda]; % xpcdata6.Tr0xF8g(102500:end) %% weight % remember q is weight between 0 and 1 (I think this will be best) - p0 is initial start value for P. %Q_wood = [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]; Q_wood = [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 10^-6 0.1*10^-3]; P0_wood = [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]; R_wood=[.1 0;0 .1];