[Scilab-users] Linearizing a control system using "lin" (porting from MATLAB)

[gabrielgga]

I am currently trying to port some of my MATLAB code to Scilab. 

To linearize odes in MATLAB, I have created a set of symbolic variables and
expressions and computed the jacobians so I could later use them in ss2tf()
and the use lsim to simulate the system's response to a series of step
inputs.

Relevant MATLAB code (the full code is at http://pastebin.ca/3727857
(comments in portuguese):

    % Linearization

    % Symbolic variables
    syms s_f_e s_T_e s_f_t s_T_t_e...
       s_Cp s_Cp_t s_rho s_rho_t s_Area s_Area_t s_V_t s_U s_Cv

    % State variables
    syms s_h s_T s_T_t

    % State arrays
    u = [s_f_e s_T_e s_f_t];
    x = [s_h s_T s_T_t];

    % States
    % h 
    F1 = (s_f_e - s_Cv*sqrt(s_h))/s_Area;
    % T
    F2 = ((s_f_e/s_h)*(s_T_e-s_T) + ...
            (s_U*s_Area_t/(s_rho*s_Cp*s_h))*(s_T_t-s_T))/s_Area;
    % T_t
    F3 = (s_f_t*(s_T_t_e-s_T_t) - ...
             (s_U*s_Area_t/(s_rho_t*s_Cp_t))*(s_T_t-s_T))/s_V_t;

    F = [F1 F2 F3];

    % Outputs
    G = [s_h s_T s_T_t];

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

    % manipulating A, B, C, D to get numeric values, using ss2tf to get
transfer functions

    % Linear system response to steps in input f_e
    dt = 0.1;
    tf = 60;
    t = 0:dt:tf;
    amp = 0.2*f_e_bar;

    % u consists of two steps at times st_f_e(1) and st_f_e(2) with
amplitude amp and -amp
    u = [zeros(st_f_e(1)/dt,1);...
    amp*ones((st_f_e(2) - st_f_e(1))/dt, 1);...
    zeros((tf - st_f_e(2))/dt + 1,1)];

    figure(1);

    % h
    y = lsim(h_fe,u,t) + h_0;

    subplot(3,1,1);
    plot(t,y,'--r');


I understand Scilab doens't support symbolic mathematics, so after a bit of
research I have this:

    // Linearization

    // State arrays

    u = [f_e_bar T_e_bar f_t_bar]';
    x = [h_0 T_0 T_t_0]';

    // States

    function [y,xdot] = f(u,x)
    
        // Process parameters
        Cp        = 0.75; // kJ/kg*K
        Cp_t      = 1;    // kJ/kg*K
        rho       = 1094; // kg/m³
        rho_t     = 1194; // kg/m³
        Area      = 2;    // m²
        Area_t    = 1.5;  // m²
        V_t       = 5;    // m³
        U         = 150;  // kJ/m²*s*K
        Cv        = 4;    // constante da válvula
        T_t_e_bar = 540;  // K
    
        // h
        xdot(1) = (u(1) - Cv*sqrt(x(1)))/Area;
        // T
        xdot(2) = ((u(1)/x(1))*(u(2)-x(2)) + ...
                (U*Area_t/(rho*Cp*x(1)))*(x(3)-x(2)))/Area;
        // T_e
        xdot(3) = (u(3)*(T_t_e_bar-x(3)) - ...
                 (U*Area_t/(rho_t*Cp_t))*(x(3)-x(2)))/V_t;
    
        y = x;
    
    endfunction

    [A,B,C,D] = lin(f,u,x)

    sys = syslin('c',A,B,C,D,x);

    // f_e
    h = ss2tf(sys);
    h_fe   = h(1,1);

    // Resposta do sistema linear (f_e)
    dt = 0.1;
    tf = 60;
    t = [0:dt:tf];
    amp = 0.2*f_e_bar;
    u = [zeros(st_f_e(1)/dt,1);...
        amp*ones((st_f_e(2) - st_f_e(1))/dt, 1);...
        zeros((tf - st_f_e(2))/dt + 1,1)]';
    
    // h
    y = csim(u,t,h_fe) + h_0;

    figure(1);
    subplot(3,1,1);
    plot(t,y,'--r');

However, the output is completely wrong. I suspect the main problem is in
the function f(u,x) I pass to lin(). Is this the way to linearize a
non-linear system of odes? Is there a simpler one?



--
View this message in context: http://mailinglists.scilab.org/Linearizing-a-control-system-using-lin-porting-from-MATLAB-tp4034786.html
Sent from the Scilab users - Mailing Lists Archives mailing list archive at Nabble.com.