Linearize non-linear system using Matlab/Simulink

188 views (last 30 days)
Wouter
Wouter on 10 Dec 2012
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.

Accepted Answer

Ryan G
Ryan G on 10 Dec 2012
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.
  2 Comments
Ryan G
Ryan G on 11 Dec 2012
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.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by