Linearize non-linear system using Matlab/Simulink

I have a non-linear system;

This system corresponds to the following blockdiagram in Simulink;

With

g       = 9.81;   % gravitational force [m/s^2]
rho     = 980;    % mass density [kg/m^3]
A       = 1;      % area [m^2]
K       = 0.01;   % valve coefficient [-]
A_v     = 0.001;  % valve cross-sectional area [m^2]
phi_o_0 = 0.001;  % initial flow out
m_0     = (phi_o_0 / (K * A_v))^2 * A / g; % initial mass [kg]

I need to linearize this system around the working point m_0. With A_v and phi_i as inputs and phi_o, m, h and p_i as outputs.

This I did easily using just plain Matlab code;

close all;
clear all;
clc;

% model variables (Area = V instead of A).
syms V rho g K;

% state variables
syms phi_i A_v h m p phi_o;

% state vectors
u = [phi_i; A_v];
x = [m]; 
y = [h; m; p; phi_o];

% non-linear system, dx(t)/dt = f(x,u,t)
F1  = rho * phi_i - rho * K * A_v * sqrt((g / V) * m);
F = [F1];

% non-linear system, y(t) = g(x,u,t)
G2 = m;
G1 = rho * V * G2;
G3 = (g * G2) / V;
G4 = K * A_v * sqrt(G3);
G = [G1; G2; G3; G4];

% compute jacobian
A.symbolic = jacobian(F, x);
B.symbolic = jacobian(F, u);
C.symbolic = jacobian(G, x);
D.symbolic = jacobian(G, u);

% Algebraic value of operating point see ex1a.
m_0 = sym(8966455680130479/8796093022208);

% compute matrices A, B, C, D
A.algebraic = simplify(subs(A.symbolic, {A_v V rho g K m}, ...
                [sym(0.01) sym(1) sym(980) sym(9.81) sym(0.01) m_0]));

B.algebraic = simplify(subs(B.symbolic, {A_v V rho g K m}, ...
                [sym(0.01) sym(1) sym(980) sym(9.81) sym(0.01) m_0]));

C.algebraic = simplify(subs(C.symbolic, {A_v V rho g K m}, ...
                [sym(0.01) sym(1) sym(980) sym(9.81) sym(0.01) m_0]));

D.algebraic = simplify(subs(D.symbolic, {A_v V rho g K m}, ...
                [sym(0.01) sym(1) sym(980) sym(9.81) sym(0.01) m_0]));             

% compute numerical values            
A.eval = eval(A.algebraic);
B.eval = eval(B.algebraic);
C.eval = eval(C.algebraic);
D.eval = eval(D.algebraic);

% linearized system
linsys = ss(A.eval, B.eval, C.eval, D.eval);

However just out of curiousity I also wanted Matlab compute the linearized system, http://www.mathworks.nl/help/slcontrol/ug/linearize.html.

I however miserably fail using the linearize functionality and I am not exactly sure what I am doing wrong. So I was hoping someone can help me to get this fixed.

Currently I have

sys = 'ex1a_model';
load_system(sys);
open_system(sys);

% Inputs
sys_io(1) = linio('ex1a_model/phi_i(t)',1,'in');
sys_io(2) = linio('ex1a_model/A_v',2,'in');

% Outputs
sys_io(3) = linio('ex1a_model/Out m(t)',1,'out');
sys_io(4) = linio('ex1a_model/Out h(t)',2,'out');
sys_io(5) = linio('ex1a_model/Out p_i(t)',3,'out');
sys_io(6) = linio('ex1a_model/Out phi_o(t)',6,'out');

% update linio
setlinio(sys,sys_io);

% Set openloop
sys_io(3).OpenLoop='on';
sys_io(4).OpenLoop='on';
sys_io(5).OpenLoop='on';
sys_io(6).OpenLoop='on';

% Linearize
linsys = linearize(sys,sys_io);

I hope someone can help me with this.

m-file: http://dl.dropbox.com/u/20782274/ex1b.m

mdl-file: http://dl.dropbox.com/u/20782274/ex1a_model.mdl

Answer 1

[Ryan G]

Your best bet is to perform the first linearization in the Linear Analysis Tool and then use the 'Generate M-File' option in the tool to create the code.

This way you know you are going to get the syntax correct and can easily change the parameters.

This is covered a bit in the trim and linearization video.

One note, you can load or open the system, you don't need to do both. One loads without opening the Simulink GUI, the other will load AND open the gui.

[Wouter]

Thanks Ryan G. Your link to the Linear Analysis Tool (LAT) helped me a lot. I didn't know there was something like that.

Funny thing is that the LAT gives me a slightly different result that what I compute using my Matlab script...

My Matlab script gives:

linsys =

    a = 
                x1
     x1  -0.004807

    b = 
           u1    u2
     x1   980  -980

    c = 
                x1
     y1    0.00102
     y2          1
     y3       9.81
     y4  4.905e-06

    d = 
         u1  u2
     y1   0   0
     y2   0   0
     y3   0   0
     y4   0   1

    Continuous-time state-space model.

But the LAT tool gives:

linsys =

    a = 
           x1
     x1     0

    b = 
           u1    u2
     x1   980  -980

    c = 
               x1
     y1   0.00102
     y2         1
     y3      9.81
     y4         0

    d = 
         u1  u2
     y1   0   0
     y2   0   0
     y3   0   0
     y4   0   1

    Continuous-time state-space model.

Where does this difference comes from? Is it due to very small differences during calculations?

[Ryan G]

Don't know for sure, there is an option you can check to 'launch diagnostic viewer'. This may help you better understand how each block is being linearized.

To determine which model you feel is more accurate, you can run a step response through the linear models and compare that to a step in the main system. If that is not enough you can try sine sweeps (available in the linear analysis tool) and compare the bode of the sine sweep analysis to the linear system analysis.