Example 4: Hydraulic System
The above figure presents the schematic diagram for a hydraulic ram. For more details of the plant and an industrial application, see also [Kugi2000] A. Kugi, Nonlinear Control Based on Physical Models, Lecture Notes in Control and Information Sciences 260, SpringerVerlag, 2000. The continuity equations of this system read as
> 
diff(P1(t),t) = Eoil (Q1  vp*A1  Cleak_int*(P1P2)Cleak_ext1*P1)/(V01+xp*A1);
diff(P2(t),t) = Eoil (Q2 + vp*A2 + Cleak_int*(P1P2)Cleak_ext2*P2)/(V02xp*A2);
diff(xp(t),t) = vp;
diff(vp(t),t) = (P1*A1  P2*A2)/m;

Here P1, P2 denote the pressures in the chambers, xp is the piston displacement measured from the center, vp is the piston velocity. Additional parameters are the cylinder volumes V01, V02 (measured at xp = 0), the piston area A, the effective bulk modulus Eoil and the valve flows Q1 and Q2. The supply pressure is given by Ps, the tank pressure by Pt, the mass of the moving parts is m. The external and internal leakage flows are assumed to be laminar with the leakage coefficients Cleak_ext1, Cleak_ext2 and Cleak_int. The valve flows for positive and negative spool displacements xs are given by the equations below, Kd is a valve parameter.
If xs > 0 the flows are:
Q1_pos := Kd*xs*sqrt(PsP1);
Q2_pos := Kd*xs*sqrt(P2Pt);
If xs < 0 the flows are
Q1_neg := Kd*xs*sqrt(P1Pt);
Q2_neg := Kd*xs*sqrt(PsP2);
If xs = 0 the flows are
Q1_zero := 0;
Q2_zero := 0;
Two different models for positive and negative valve displacements are necessary, the plant input is the valve position xs and the plant output is the piston position xp.
> 
f:=[Eoil/(V01+xp*A1)*(vp*A1 Cleak_int*(P1P2)Cleak_ext1*P1),Eoil/(V02xp*A2)*(vp*A2+ Cleak_int*(P1P2)Cleak_ext2*P1),vp,(P1*A1P2*A2)];
g_pos:=[Eoil/(V01+xp*A1)*Kd*sqrt(PsP1),Eoil/(V02xp*A2)*Kd*sqrt(P2Pt),0,0];
g_neg:=[Eoil/(V01+xp*A1)*Kd*sqrt(P1Pt),Eoil/(V02xp*A2)*Kd*sqrt(PsP2),0,0];
xx:=[P1,P2,xp,vp];

Often in hydraulic systems the assumption PsP1 = P2  Pt is made, in this case the flow Q = Q1 = Q2 can be used as input and the vector field g becomes
> 
g_Q:=[Eoil/(V01+xp*A1),Eoil/(V02xp*A2),0,0];

> 
`AIsys/Global/simp` := proc(x)
simplify(x);
end:

Setup the mathematical model for the AIsys package
> 
f:=M_AIsys(f);
g_pos:=M_AIsys(g_pos);
g_neg:=M_AIsys(g_neg);
g_Q:=M_AIsys(g_Q);

Computation of the reachability distribution, to check if the systems hangs up we want more outputs
> 
print(`For positive valve displacement`);
RDistr(f,g_pos,xx):

`Initialization: delta:={}`
`delta := delta + g1`
`Computing Lie bracket: [f,g1]`
`delta := delta + [f,g1]`
`Computing Lie bracket: [delta2,delta1]`
`delta := delta + [delta2,delta1]`
`Computing Lie bracket: [f,delta2]`
`delta := delta + [f,delta2]`
`Dimension of reachability distribution:`, 4
More outputs can be obtained by setting infolevel[AIsys]:=3: or infolevel[AIsys]:=4:
What happens if the flow Q is used as input?
`Initialization: delta:={}`
`delta := delta + g1`
`Computing Lie bracket: [f,g1]`
`delta := delta + [f,g1]`
`Computing Lie bracket: [delta2,delta1]`
`Computed Lie bracket is linear dependent.`
`Computing Lie bracket: [f,delta2]`
`delta := delta + [f,delta2]`
`Computing Lie bracket: [delta3,delta1]`
`Computed Lie bracket is linear dependent.`
`Computing Lie bracket: [delta3,delta2]`
`delta := delta + [delta3,delta2]`
`Dimension of reachability distribution:`, 4
Is the system still reachable without any leakages? To answer this question the siderelations in the simplifiers routines are adjusted in the following way:
> 
`AIsys/SIMP`[2]:={Cleak_ext1=0,Cleak_ext2=0,Cleak_int=0};

`Initialization: delta:={}`
`delta := delta + g1`
`Computing Lie bracket: [f,g1]`
`delta := delta + [f,g1]`
`Computing Lie bracket: [delta2,delta1]`
`Computed Lie bracket is linear dependent.`
`Computing Lie bracket: [f,delta2]`
`delta := delta + [f,delta2]`
`Computing Lie bracket: [delta3,delta1]`
`Computed Lie bracket is linear dependent.`
`Computing Lie bracket: [delta3,delta2]`
`Computed Lie bracket is linear dependent.`
`Computing Lie bracket: [f,delta3]`
`Computed Lie bracket is linear dependent.`
`Dimension of reachability distribution:`, 3
This result coincides with the analysis of the linearized system.
Next the observability is checked, readjusting the simplifier routines.
> 
print(`For negative valve displacement`);
ODistr(f,g_neg,xp,xx):

`Initialization: delta := {}`
`Computing differential dh`
`delta := delta + dh`
`Computing Lie derivative: [delta1,f]`
`delta := delta + Lie derivative`
`Computing Lie derivative: [delta1,g1]`
`Computed Lie derivative is linear dependent.`
`Computing Lie derivative: [delta2,f]`
`delta := delta + Lie derivative`
`Computing Lie derivative: [delta2,g1]`
`Computed Lie derivative is linear dependent.`
`Computing Lie derivative: [delta3,f]`
`delta := delta + Lie derivative`
`Dimension of observability codistribution:`, 4
What happens if Q is used as input?
> 
print(`With flow as input`);
ODistr(f,g_Q,xp,xx):

`Initialization: delta := {}`
`Computing differential dh`
`delta := delta + dh`
`Computing Lie derivative: [delta1,f]`
`delta := delta + Lie derivative`
`Computing Lie derivative: [delta1,g1]`
`Computed Lie derivative is linear dependent.`
`Computing Lie derivative: [delta2,f]`
`delta := delta + Lie derivative`
`Computing Lie derivative: [delta2,g1]`
`Computed Lie derivative is linear dependent.`
`Computing Lie derivative: [delta3,f]`
`delta := delta + Lie derivative`
`Dimension of observability codistribution:`, 4
What happens without any leakages?
> 
`AIsys/SIMP`[2]:={Cleak_ext1=0,Cleak_ext2=0,Cleak_int=0}:

> 
print(`With flow as input and no leakages`);
ODistr(f,g_Q,xp,xx):

`Initialization: delta := {}`
`Computing differential dh`
`delta := delta + dh`
`Computing Lie derivative: [delta1,f]`
`delta := delta + Lie derivative`
`Computing Lie derivative: [delta1,g1]`
`Computed Lie derivative is linear dependent.`
`Computing Lie derivative: [delta2,f]`
`delta := delta + Lie derivative`
`Computing Lie derivative: [delta2,g1]`
`Computed Lie derivative is linear dependent.`
`Computing Lie derivative: [delta3,f]`
`Computed Lie derivative is linear dependent.`
`Computing Lie derivative: [delta3,g1]`
`Computed Lie derivative is linear dependent.`
`Dimension of observability codistribution:`, 3
Next some other properties are checked
Computing relative degree with MIMO and SISO algorithms
> 
MRelDeg(f,[g_pos],[xp],xx):
print(`SISO Algorithm`);
RelDeg(f,g_pos,xp,xx,true):

`Computing Lie derivative: L[g1]h1`
`Computed Lie derivative is 0`
`Computing Lie derivative: lie := L[f]^1h1`
`Computing Lie derivative: lie := L[g1]lie `
`Computed Lie derivative is 0`
`Computing Lie derivative: lie:=L[f]^2h1`
`Computing Lie derivative: lie := L[g1]lie `
`Lie derivative <> 0 > r1 :=3`
`Checking the rank of the decoupling matrix`
`Decoupling matrix seems to be regular.`
`The relative degree vector is `, [3]
`Computing Lie derivative : L[g]L[f]h`
`Computing Lie derivative of output h along f: L[f]^2h`
`Checking column dimension of jacobian A = d(t_vec)/dx`
`Computed relative degree is 3`, ` jacobian A is regular`
Computing the maximum relative degree
> 
MaxRelDeg(f,[g_pos],xx);

`Initialization: delta := {}`
`delta := delta + g1`
`Computing involutive closure ...`
`Rank of involutive closure of delta is 1`
`Compute Lie bracket: [f,g1]`
`delta := delta + [f,g1]`
`Computing involutive closure ...`
`Rank of involutive closure of delta is 3`
`Computing Lie bracket [f,delta2]`
`delta := delta + [f,delta2]`
`Computing involutive closure ...`
`Rank of involutive closure of delta is 4`
So the maximum relative degree is 3
For control purposes the nonlinear feedback for the Input/Output Linearization can be computed by UTrans
If xs > 0
print(`For positive valve displacement use transformation:`);
nlp:=UTrans(f,g_pos,xp,xx,v,true);
If xs < 0
print(`For negative valve displacement use transformation:`);
nln:=UTrans(f,g_neg,xp,xx,v,true);
`Computing Lie derivative : L[g]L[f]h`
`Computing Lie derivative of output h along f: L[f]^2h`
`Checking column dimension of jacobian A = dt/dx`
`Computed relative degree is 3 and jacobian of A is regular`
`Used coefficients: a0, a1, a2, `
`Coefficient: a3 set to 1`
`Computing Lie derivative : L[g]L[f]h`
`Computing Lie derivative of output h along f: L[f]^2h`
`Checking column dimension of jacobian A = dt/dx`
`Computed relative degree is 3 and jacobian of A is regular`
`Used coefficients: a0, a1, a2, `
`Coefficient: a3 set to 1`
The problem in here is that the information of the whole state is necessary to implement this feedback law. Usually, only the piston position xp and the pressures P1 and P2 are measured. The velocity signal vp can only be obtained by approximate differentiation of xp. But this procedure often fails, since in many industrial environments the transducer noise is considerable high. Therefore, we pose the question if it is possible to find an exact linearizing output such that the static feedback law is independent of the velocity vp? This question can be answered by the algorithm
ElimSt
, see, e.g., [Schlacher2001] K. Schlacher, A. Kugi, R. Novak, Input to Output Linearization with Constrained Measurements, In Preprints of the 5th IFAC Symposium on "Nonlinear Control Systems", Saint Petersburg, Russia, July 46, 2001, pp.471476, 2001.
To get less information we set
The measurements are
The computation of the restrictions gives
> 
pdes:=ElimSt(f,g_pos,xx,meas,h);

> 
pde1:=pdes[1]:
pde2:=pdes[2]:

The solution of the first restriction is clear
> 
restr1:=pdsolve(pde1,h(op(xx)));

Substituting this result in the second one gives the desired output
> 
restr2:=simplify(subs(pde1=0,h(op(xx))=h(op(meas)),pde2));

> 
output:=pdsolve(restr2,h(op(meas)));

For instance the output
> 
z:=P1*A1P2*A2+Eoil*ln((V01+A1*xp)^A1/(V02A2*xp)^A2);

gives the desired result.
> 
simplify(subs(h(op(xx))=z,pde1));
simplify(subs(h(op(xx))=z,pde2));

The nonlinear feedback for positive valve displacement is now independent of the velocity and becomes
> 
UTrans(f,g_pos,z,xx,v);

Sometimes only the load pressure P1  P2 can be measured, what happens in this case?
> 
pdes:=ElimSt(f,g_pos,xx,meas,h);

"Problem has no solution."
Thus, in the general case the problem has no solution. But what happens without leakages and with the flow Q as input?
> 
`AIsys/SIMP`[2]:={Cleak_ext1=0,Cleak_ext2=0,Cleak_int=0}:

> 
pdes:=ElimSt(f,g_Q,xx,meas,h);

> 
pde1:=pdes[1]:
pde2:=pdes[2]:
pde3:=pdes[3]:

The influence of the first and the second restriction is clear
> 
pdsolve(pde1,h(op(xx)));

> 
pdsolve(pde2,h(op(xx)));

Thus, we get the same output as before
> 
pdsolve(subs(h(op(xx))=h(P1,P2,xp),pde3),h(P1,P2,xp));
