Solvin Differential Algebraic Equations

Please how can I solve a DAE for a gas lift system network. I have 5 differential equations, and 18 algebraic equations. I tried following the codes to create mine but it is not running. Here is a sample:

function out = gas_lift_ftn(t,y)

%Define other parameters%

mu_o = 0.01;

pi = 3.14;

Lw = 1500;

Hw = 1000;

Dw = 0.121;

Lbh = 500;

Hbh = Lbh;

Dbh = Dw;

La = Lw;

Ha = Hw;

Da = 0.189;

Lr = Lbh;

Hr = Lbh;

Dr = Dw;

Aw = (pi/4)*Dw^2;

Abh = Aw;

Ar = Aw;

Aa = (pi/4)*Da;

Va = La*pi/4*(Da^2-Dw^2);

rho_o = 800;

Civ = 0.1e-03; Crh = 10e-3; Cpc = 2e-3;

Pr = 150; Ps = 20; PI = 0.7;

Ta = 28+273; Tw = 32+273; Tr = 30+273; Mw = 20; R = 0.083145;

g = 9.8*10e-5; GOR = 0.1; wgl = 2.217;

out = [ wgl - y(15)

y(15) - y(19) + y(18)

y(20) - y(21)

y(19) - y(22)

y(21) - y(23)

((((Ta*R)/(Va*Mw))+((g*La)/(La*Aa)))*y(1)) - y(6)

(Mw*y(6)) -(Ta*R*y(9))

(((((Tw*R)/Mw))*((y(2))/((Lw*Aw)+(Lbh*Abh)-(y(3)/rho_o)))- 0.5*((y(2)+(y(3))/(Lw*Aw))*g*Hw))-y(7))

y(2)+ y(3)-(rho_o*Lbh*Abh)-(Lw*Aw*y(10))

(y(7)+((g/(Lw*Aw))*(y(3)+y(2)-(rho_o*Lbh*Abh)*Hw))+((128*mu_o*Lw*y(16))/(pi*Dw^4*((y(2) + y(3))*y(7)*Mw*rho_o)/(y(3)*y(7)*Mw + rho_o*R*Tw*y(2)))))-y(8)

(y(8)+ (y(10)*g*Hbh)+((128*mu_o*Lbh*y(20))/(pi*Dbh^4*rho_o)))-y(12)

(Tr*R*y(4))-(Mw*Lr*Ar*y(13))

y(4)+ y(5)-(Lr*Ar*y(11))

(y(13) + y(11)*g*Hr + ((128*mu_o*Lr*y(17))/(pi*Dr^4*((y(4) + y(5))*y(13)*Mw*y(11))/(y(5)*y(13)*Mw + y(11)*R*Tr*y(4)))) - y(14))

(Civ*sqrt(y(9)*(y(6)-y(8))))-y(15)

(Cpc*sqrt(y(10)*(y(7)-y(14))))-y(16)

(Crh*sqrt(y(11)*(y(6)-Ps)))- y(17)

(PI*(Pr -y(12)))- y(20)

GOR* y(20) -y(18)

(y(2)*y(16)) - (y(2)*y(19)) - (y(3)*y(19))

(y(3)*y(16)) - (y(2)*y(21)) - (y(3)*y(21))

(y(4)*y(17)) - (y(4)*y(22)) - (y(4)*y(22))

(y(5)*y(17)) - (y(5)*y(23)) + (y(5)*y(23))];

end

%gaslift systems network

tspan = [0:1:60];

% The script for gas lift network model defining initial conditions and

% boundaries

y0 =[2812.3; 1269.75; 9200.4; 136.9; 2300.2; 160; 80.52; 113.87; 127.9;

340.29; 428.8; 130.54; 30; 50.77; 1.853;63.65; 7.71; 55.93; 13.62; 1.36;

205.86; 11.56; 194.3];

% call

[t,y] = ode15s(@gas_lift_ftn,tspan,y0);

disp('y(1), y(2), y(3), y(4), y(5)')

disp(y(5:1))

[Torsten]

And where do you define which equations are differential equations and which are algebraic equations (i.e. the mass matrix M for your DAE system) ?

[Patience Shamaki]

Actually that is my challenge assigning the variables. I am not familiar with solving DAE's and most I see requires me to differentiate the algebraic equations into full ODE's or just a few etc. I just need guidance on how to go about the coding.

[Torsten]

Then please let me know which of your 23 equations read "0 = out(i)" and which read "dy(i)/dt = out(i)".

[Patience Shamaki]

Differential equations: (dx(1)/dt = wgl - wiv; dx(2)/dt = wiv - wpg+ wrg .....)

dy(1)/dt= wgl - y(15);

dy(2)/dt= y(15) - y(19) + y(18);

dy(3)/dt= y(20) - y(21);

dy(4)/dt= y(19) - y(22);

dy(5)/dt= y(21) - y(23);

0= ((((Ta*R)/(Va*Mw))+((g*La)/(La*Aa)))*y(1)) - y(6);

0= (Mw*y(6)) -(Ta*R*y(9));

0=(((((Tw*R)/Mw))*((y(2))/((Lw*Aw)+(Lbh*Abh)-(y(3)/rho_o)))- 0.5*((y(2)+(y(3))/(Lw*Aw))*g*Hw))-y(7));

0= y(2)+ y(3)-(rho_o*Lbh*Abh)-(Lw*Aw*y(10));

0= (y(7)+((g/(Lw*Aw))*(y(3)+y(2)-(rho_o*Lbh*Abh)*Hw))+((128*mu_o*Lw*y(16))/(pi*Dw^4*((y(2) + y(3))*y(7)*Mw*rho_o)/(y(3)*y(7)*Mw + rho_o*R*Tw*y(2)))))-y(8);

0= (y(8)+ (y(10)*g*Hbh)+((128*mu_o*Lbh*y(20))/(pi*Dbh^4*rho_o)))-y(12);

0= (Tr*R*y(4))-(Mw*Lr*Ar*y(13));

0= y(4)+ y(5)-(Lr*Ar*y(11));

0= (y(13) + y(11)*g*Hr + ((128*mu_o*Lr*y(17))/(pi*Dr^4*((y(4) + y(5))*y(13)*Mw*y(11))/(y(5)*y(13)*Mw + y(11)*R*Tr*y(4)))) - y(14));

0= (Civ*sqrt(y(9)*(y(6)-y(8))))-y(15);

0= (Cpc*sqrt(y(10)*(y(7)-y(14))))-y(16);

0= (Crh*sqrt(y(11)*(y(6)-Ps)))- y(17);

0= (PI*(Pr -y(12)))- y(20);

0= GOR* y(20) -y(18);

0= (y(2)*y(16)) - (y(2)*y(19)) - (y(3)*y(19));

0= (y(3)*y(16)) - (y(2)*y(21)) - (y(3)*y(21));

0= (y(4)*y(17)) - (y(4)*y(22)) - (y(4)*y(22));

0=(y(5)*y(17)) - (y(5)*y(23)) + (y(5)*y(23));

Answer 1

[Torsten]

function main
  %gaslift systems network
  tspan = [0:1:60];
  % The script for gas lift network model defining initial conditions and
  % boundaries
  y0 =[2812.3; 1269.75; 9200.4; 136.9; 2300.2; 160; 80.52; 113.87; 127.9; 
      340.29; 428.8; 130.54; 30; 50.77; 1.853;63.65; 7.71; 55.93; 13.62; 1.36;
      205.86; 11.56; 194.3];
  M = diag([ones(1,5) zeros(1,18)]);
  options = odeset('Mass',M);
  % call 
  [t,y] = ode15s(@gas_lift_ftn,tspan,y0,options);
end

function out = gas_lift_ftn(t,y)
  %Define other parameters%
  mu_o = 0.01; 
  pi = 3.14; 
  Lw = 1500;
  Hw = 1000; 
  Dw = 0.121;
  Lbh = 500; 
  Hbh = Lbh; 
  Dbh = Dw; 
  La = Lw; 
  Ha = Hw; 
  Da = 0.189;
  Lr = Lbh;
  Hr = Lbh; 
  Dr = Dw; 
  Aw = (pi/4)*Dw^2;
  Abh = Aw;
  Ar = Aw;
  Aa = (pi/4)*Da;
  Va = La*pi/4*(Da^2-Dw^2);
  rho_o = 800;
  Civ = 0.1e-03; Crh = 10e-3; Cpc = 2e-3;
  Pr = 150; Ps = 20; PI = 0.7; 
  Ta = 28+273; Tw = 32+273; Tr = 30+273; Mw = 20; R = 0.083145;
  g = 9.8*10e-5; GOR = 0.1; wgl = 2.217;

  out=zeros(23,1);
  out(1) = wgl - y(15);
  out(2) = y(15) - y(19) + y(18);
  out(3) = y(20) - y(21);
  out(4) = y(19) - y(22);
  out(5) = y(21) - y(23);
  out(6) = ((((Ta*R)/(Va*Mw))+((g*La)/(La*Aa)))*y(1)) - y(6);
  out(7) = (Mw*y(6)) -(Ta*R*y(9));
  out(8) = (((((Tw*R)/Mw))*((y(2))/((Lw*Aw)+(Lbh*Abh)-(y(3)/rho_o)))- 0.5*((y(2)+(y(3))/(Lw*Aw))*g*Hw))-y(7));
  out(9) = y(2)+ y(3)-(rho_o*Lbh*Abh)-(Lw*Aw*y(10));
  out(10) = (y(7)+((g/(Lw*Aw))*(y(3)+y(2)-(rho_o*Lbh*Abh)*Hw))+((128*mu_o*Lw*y(16))/(pi*Dw^4*((y(2) + y(3))*y(7)*Mw*rho_o)/(y(3)*y(7)*Mw + rho_o*R*Tw*y(2)))))-y(8);
  out(11) = (y(8)+ (y(10)*g*Hbh)+((128*mu_o*Lbh*y(20))/(pi*Dbh^4*rho_o)))-y(12);
  out(12) = (Tr*R*y(4))-(Mw*Lr*Ar*y(13));
  out(13) = y(4)+ y(5)-(Lr*Ar*y(11));
  out(14) = (y(13) + y(11)*g*Hr + ((128*mu_o*Lr*y(17))/(pi*Dr^4*((y(4) + y(5))*y(13)*Mw*y(11))/(y(5)*y(13)*Mw + y(11)*R*Tr*y(4)))) - y(14));
  out(15) = (Civ*sqrt(y(9)*(y(6)-y(8))))-y(15);
  out(16) = (Cpc*sqrt(y(10)*(y(7)-y(14))))-y(16);
  out(17) = (Crh*sqrt(y(11)*(y(6)-Ps)))- y(17);
  out(18) = (PI*(Pr -y(12)))- y(20);
  out(19) = GOR* y(20) -y(18);
  out(20) = (y(2)*y(16)) - (y(2)*y(19)) - (y(3)*y(19));
  out(21) = (y(3)*y(16)) - (y(2)*y(21)) - (y(3)*y(21))
  out(22) = (y(4)*y(17)) - (y(4)*y(22)) - (y(4)*y(22))
  out(23) = (y(5)*y(17)) - (y(5)*y(23)) + (y(5)*y(23))];
end

[Patience Shamaki]

Thanks this was very helpful. I am still getting this error though

Error using daeic12 (line 76)

This DAE appears to be of index greater than 1.

Error in ode23t (line 288)

[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in gaslift_script (line 11)

[t,y] = ode15s(@gas_lift_ftn,tspan,y0,options);

[Torsten]

Then nothing helps but rechecking whether you correctly set up the DAE system.

Explicity write down which variable you want to solve with which equation,e.g.

y(1) with equation 1

y(2) with equation 2

y(3) with equation 3

y(4) with equation 4

y(5) with equation 5

y(6) with equation 6

y(9) with equation 7

y(7) with equation 8

....

Do you get an equation for every unknown ?

[Patience Shamaki]

Yes I have an equation for every unknow but they are mostly dependent on the previous or the next, I got my initial conditions of the algebraic eqautions by making some assumptions based the process.

[Jaffar Ali Lone]

@Torsten How to solve discrete time descriptor system of the following type

Ex(k) =Ax(k) + Bu(k) where E is Singular matrix

For example

E = [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 0 0;0 0 0 0 0 0];

A = [0 0 0 0 0.43478 0;0 -0.16666 0 0 0.000666 0;0 0 0 0 0 0.5;0 0 0 -0.14285 0 0.0005;22.9676 -21.9676 0 -1 0.0025 -0.0015;0 0 0 0 1 1];

C = [22.9676 1 0 0 0.0025 0;0 22.9676 0 1 0 0.0015];

B = [0;0;0;0;1;1];

u = sin(t)

[Torsten]

x{k}=(E-A)\(B*u(k))