No answers, because you've just taught me about the existence of 'lin'
-- but here are some thoughts and questions:

First - I don't see where you're putting parameter values into Matlab's
symbolic variables to get an answer.  Is this done interactively or
something?

Second - I just tried lin with a very simple system (xdot = (x + u)^2, y
= x).  It seems to work, although it didn't get the values exactly and I
don't see a way to control the precision of the numerical derivative
that gets taken under the hood.

It may be productive to try out some simpler systems and see where
you're going wrong.  Alternately, to test your function definition, you
could modify it to work with the ODE solver and see if you get the
expected behavior out of it.

Third - What's "completely wrong" -- what should the thing be doing, and
what is it doing?  I assume that Matlab is working correctly, with some
cross-check to reality, and Scilab isn't.

If you don't mind working with two applications, Maxima would do the
symbolic calculations that you're calling out, and the answers could be
transferred to Scilab.  Unfortunately the syntax is considerably
different, so you'd have two learning curves instead of one.  At one
point there was an integration of Maxima into Scilab, but I don't know
how well it worked, or if it's still being supported.

I haven't gone this far down the automation road with Scilab, but I'm
99.44% sure that you can write a Maxima script that does the symbolic
calculation and spits out some Scilab code, then call Maxima from
Scilab, execute the Scilab code, and have a system.  It's much more
horribly indirect than what Matlab appears to be doing for you, though.

On Wed, 2016-10-12 at 10:00 -0700, gabrielgga wrote:
> 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.
> _______________________________________________
> users mailing list
> users@lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users

-- 

Tim Wescott
www.wescottdesign.com
Control & Communications systems, circuit & software design.
Phone: 503.631.7815
Cell:  503.349.8432


_______________________________________________
users mailing list
users@lists.scilab.org
http://lists.scilab.org/mailman/listinfo/users

Reply via email to