Hello everyone, thanks so much for all the help I am getting here to learn MATLAB. Please how can I convert my script into a function. I did an integration using Euler explicit method and I want to use it for the objective function of my optimization, so I want to convert the script into a function. Please I really need help with this.
%integration of gas-lift ODEs using euler explicit methods
% clear all
% clc
% close
%timesteps
tf = 3600;
h = 0.1;
t = 0:h:tf;
n = length(t);
%initialising
x10 = 1.0568;
x20 = 1.0606;
x30 = 0.6000;
x40 = 0.6700;
x50 = 6.0000;
x60 = 5.5047;
x70 = 0.1265;
x80 = 0.9863;
x(1,1) = x10;
x(2,1) = x20;
x(3,1) = x30;
x(4,1) = x40;
x(5,1) = x50;
x(6,1) = x60;
x(7,1) = x70;
x(8,1) = x80;
%parameters
L_w = 1500;
H_w = 1000;
D_w = 0.121;
L_bh = 500;
H_bh = 500;
D_bh = 0.121;
L_a = L_w;
H_a = H_w;
D_a = 0.189;
V_a = L_a*pi*((D_a/2)^2 - (D_w/2)^2);
A_w = pi*(D_w/2)^2;%[m2]
A_bh = pi*(D_bh/2)^2;%[m2]
L_r = 500;
H_r = 500;
D_r = 0.121;
A_r = pi*(D_r/2)^2;
rho_o = 8*1e2;
rho_ro = 8*1e2;
C_iv = 0.1e-3;
C_pc = 2e-3;
C_rh = 10e-3;
GOR1 = 0.1;
GOR2 = 0.12;
dp1 = 1e-6;
dp2 = 1e-6;
p_res1 = 150;
p_res2 = 155;
PI = 0.7;
T_a = 28+273;
T_w = 32+273;
T_r = 30+273;
p_s = 20;
% p_m = p_s;
mu_o = 0.001;
slip = 1;
R = 8.3145;
Mw = 0.020;
g = 9.8;
w_gl1 =1;
w_gl2 =1;
H_k = [0 0 0 0 0 0 20642.89 -374.645 0 0; 0 0 0 0 0 0 3674.553 -163.461 0 0];
w_k = [0.001 0.001 0.001 0.001 0.02 0.02 0.01 0.03 0.001 0.001]' ;
v_k = rand(2,1);
q = cov(w_k);
r = cov(v_k);
R_k = [0.01 0; 0 0.01];
lb = [10*1e-5;10*1e-5;10*1e-5;10*1e-5;10*1e-5;10*1e-5;10*1e-5;10*1e-5;10*1e-2;10*1e-2];
ub = [10e7;10e7;10e7;10e7;10e7;10e7;10e7;10e7;50e4;50e4];
QgMax = 10;
qGLMax = 40;
a = 1;
b = 0.5;
A = [];
b = [];
Aeq = [];
beq = [];
for i=1:n
%auxiliary algebraic equations
p_a1 = ((R*T_a)/(V_a*Mw) + (g*H_a/V_a))*x10*1e3;%pressure of gas in the annulus
p_a2 = ((R*T_a)/(V_a*Mw) + (g*H_a/V_a))*x20*1e3;
rho_a1 = (Mw/(R*T_a))*p_a1;%density of gas in the annulus
rho_a2 = (Mw/(R*T_a))*p_a2;
p_wh1 = ((R*T_w/Mw)*(x30*1e3/(L_w*A_w + L_bh*A_bh - x50*1e3/rho_o))) - (((x30*1e3+x50*1e3)/(L_w*A_w))*g*H_w/2);%pressure at the wellhead
p_wh2 = ((R*T_w/Mw)*(x40*1e3/(L_w*A_w + L_bh*A_bh - x60*1e3/rho_o))) - (((x40*1e3+x60*1e3)/(L_w*A_w))*g*H_w/2);
rho_m1 = (((x30*1e3+x50*1e3)*p_wh1*Mw*rho_o)/(x50*1e3*p_wh1*Mw + rho_o*R*T_w*x30*1e3));
rho_m2 = (((x40*1e3+x60*1e3)*p_wh2*Mw*rho_o)/(x60*1e3*p_wh2*Mw + rho_o*R*T_w*x40*1e3));
p_rh = (R*T_r/Mw)*(x70*1e3/(L_r*A_r)) - (((x70*1e3+x80*1e3)/(L_r*A_r))*g*H_r/2); %25 Pressure at the riser head 1[pa]
rho_r = ((x70*1e3+x80*1e3)*p_rh*Mw*rho_ro)/(x80*1e3*p_rh*Mw + rho_ro*R*T_r*x70*1e3); % density of mixture at the riser [pa]
w_pr = C_rh*sqrt(rho_r*(p_rh - p_s*1e5)); %28 total mass flow rate through the production choke1[kg/s]
p_m = p_rh + (g/(A_r*L_r)*(x80*1e3+x70*1e3)*H_r) + (128*mu_o*L_r*w_pr/(pi*D_r^4*((x70*1e3 + x80*1e3)*p_rh*Mw*rho_ro)/(x80*1e3*p_rh*Mw + rho_ro*R*T_r*x70*1e3)));%27%Pressure at the manifold[pa]
w_pc1 = C_pc*sqrt(rho_m1*(max(1e-3,(p_wh1 - p_m))));
w_pc2 = C_pc*sqrt(rho_m2*(max(1e-3,(p_wh2 - p_m))));
w_pg1 = (x30*1e3/(max(1e-3,(x30*1e3+x50*1e3))))*w_pc1;
w_pg2 = (x40*1e3/(max(1e-3,(x40*1e3+x60*1e3))))*w_pc2;
w_po1 = (x50*1e3/(max(1e-3,(x30*1e3+x50*1e3))))*w_pc1;
w_po2 = (x60*1e3/(max(1e-3,(x40*1e3+x60*1e3))))*w_pc2;
p_w1 = p_wh1 +((g/(A_w*L_w))*(max(1e-3,(x30*1e3+x50*1e3-rho_o*L_bh*A_bh)))*H_w)+ (128*mu_o*L_w*w_pc1/(pi*D_w^4*(x30*1e3+x50*1e3)*p_wh1*Mw*rho_o)/(x50*1e3*p_wh1*Mw + rho_o*R*T_w*x30*1e3)); %((m_gt1 + m_ot1)*p_wh1*Mw.*rho_o)./(m_ot.*1e3.*p_wh.*1e5.*Mw + rho_o.*R.*T_w.*m_gt.*1e3)));;
p_w2 = p_wh2 +((g/(A_w*L_w))*(max(1e-3,(x40*1e3+x60*1e3-rho_o*L_bh*A_bh)))*H_w)+ (128*mu_o*L_w*w_pc2/(pi*D_w^4*(x40*1e3+x60*1e3)*p_wh2*Mw*rho_o)/(x60*1e3*p_wh2*Mw + rho_o*R*T_w*x40*1e3));
w_iv1 = C_iv*sqrt(rho_a1*(max(1e-3,(p_a1 - p_w1))));
w_iv2 = C_iv*sqrt(rho_a2*(max(1e-3,(p_a2 - p_w2))));
p_bh1 = p_w1 + (rho_o*g*H_bh) + (128*mu_o*L_bh*w_po1/(pi*D_bh^4*rho_o));
p_bh2 = p_w2 + (rho_o*g*H_bh) + (128*mu_o*L_bh*w_po2/(pi*D_bh^4*rho_o));
w_ro1 = PI*1e-5*(p_res1*1e5 - p_bh1);
w_ro2 = PI*1e-5*(p_res2*1e5 - p_bh2);
w_rg1 = GOR1*w_ro1;
w_rg2 = GOR2*w_ro2;
w_to = (x80*1e3/(x70*1e3+x80*1e3))*w_pr; % 29 mass flow rate of oil from the reservoir2 [kg/s]
w_tg = (x70*1e3/(x70*1e3+x80*1e3))*w_pr; %30 mass flow rate of gas from the reservoir [kg/s]
x(1,i+1) = x(1,i) + h*(w_gl1 - w_iv1)*1e-3;
x(2,i+1) = x(2,i) + h*(w_gl2 - w_iv2)*1e-3;
x(3,i+1) = x(3,i) + h*(w_iv1 + w_rg1 - w_pg1)*1e-3;
x(4,i+1) = x(4,i) + h*(w_iv2 + w_rg2 - w_pg2)*1e-3;
x(5,i+1) = x(5,i) + h*(w_ro1 - w_po1)*1e-3;
x(6,i+1) = x(6,i) + h*(w_ro2 - w_po2)*1e-3;
x(7,i+1) = x(7,i) + h*((w_pg1+w_pg2)-w_tg)*1e-3;
x(8,i+1) = x(8,i) + h*((w_po1+w_po2)-w_to)*1e-3;
x_hat_aug = [x(1,i+1);x(2,i+1);x(3,i+1);x(4,i+1);x(5,i+1);x(6,i+1);x(7,i+1);x(8,i+1);GOR1;GOR2];
h_k = (H_k*x_hat_aug);
P_0 = diag(ones(size(x_hat_aug)));
y_meas = h_k + 1e-3*h_k*randn;
P_k =F_k*P_0*F_k'+Q;
K_k = P_k*H_k'*inv(H_k*P_k*H_k'+R_k);
x_bark = x_hat_aug + K_k*(y_meas - h_k);
l=eye(10);
P_kk = (l-K_k*H_k)*P_k;
x_est =x_bark;
X_est_all(i,:)=x_est;
y_est_all(i,:)=y_meas;
h_k_all(i,:)= h_k;
% % y1 = w_to(:,end);
% y2 = w_tg(:,end);
p = [x_est(9) x_est(10)];
all_p(i,:)= p;
x_aug = [x(1,i);x(2,i);x(3,i);x(4,i);x(5,i);x(6,i);x(7,i);x(8,i);GOR1;GOR2];
x10 = x(1,i+1);
x20 = x(2,i+1);
x30 = x(3,i+1);
x40 = x(4,i+1);
x50 = x(5,i+1);
x60 = x(6,i+1);
x70 = x(7,i+1);
x80 = x(8,i+1);
