Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

Consider the following system \begin{align} \dot{x}(t)&=\sum_{i=1}^{2}\rho_{i}(x(t))\left[A_{i}x(t)+B_{i}u(t)\right]\\ y(t)&=Cx(t) \end{align}

with:

a1 = {{-3, 2}, {-0.25, 1}}
a2 = {{-1.9, -0.4}, {-2.24, -4.7}}

b1 = {{0.25}, {1}}
b2 = {{-2.5}, {1}}

c = {{1, 0.5}, {0, 1}}

where

$$\rho_1(x_1(t))=\frac{1-\tanh(x_1(t))}{2},\quad\rho_2(x_1(t))=1-\rho_1(x_1(t))$$

How do I build the state-space LPV system? The response, for example, to a sinusoidal input?

I want to translate this code from MATLAB to Mathematica

A(:,:,1)=[-3 2; -0.25 1];
A(:,:,2)=[-1.9 -0.4 ; -2.24 -4.7];
B(:,:,1)=[0.25;1] ; B(:,:,2)=[-0.25; 1];
C=[1 0.05; 0 1];%
xx=[0.1; 0]; %initial conditions
mus(1,:)=[0.5 0.5 1]; %initial rhos

%%simulation loop

for k=1:tmax/Te

    %% System    
    t(k+1)=t(k)+Te; %%time vector;

    mus(1)=(1-tanh(xx(1,k)))/10  ;  % weigting functions
    mus(2)=1-mus(1);

    u(k+1)=10*sin(t(k)); %%input
    Aa=mus(1)*A(:,:,1)+mus(2)*A(:,:,2);
    Ba=mus(1)*B(:,:,1)+mus(2)*B(:,:,2);

    %%Euler to solve the ODES
    xx(:,k+1)=xx(:,k)+Te*(Aa*xx(:,k)+Ba*u(k) );
    y(:,k+1)=C*xx(:,k+1);  %output
end

plot(t,mus,t,mus(:,1)+mus(:,2) );
plot(t,y); 
share|improve this question
    
What points (in $x$ and $u$) do you wish to linearize about? – bill s Jun 14 '13 at 9:02
    
What's $x_1$? And why is the argument of $\rho_i$ given by $x(t)$ in the first equation and by $t$ in the last equations? – sebhofer Jun 14 '13 at 9:02
    
Do you mean a step function by $u(t)$? If you do, I think it's best to be explicit! – Vincent Tjeng Jun 14 '13 at 10:51
    
$x(t)=[x_1(t)\, x_2(t)]^T$ are the states, $u$ is the input, $y(t)$ the output, and $\rho_(i)$ are the gain scheduling functions depending on the state $x_1(t)$. The system is a linear system, in fact is linear state-space parameter varying system. The idea is to simulate the response to an input $u(t)$ (for example a sinusoidal or step) similar to an standard LTI (linear time invariant system). I Matlab I made this by solving the LTI system in a for loop. But in Mathematica I don't known how to do it?. $A,B,C$ are constant matrices, $\rho_i$ changes for each step time. – user70012 Jun 14 '13 at 13:09
    
This is the code for matlab – user70012 Jun 14 '13 at 13:14

StateSpaceModel will linearize the equations. The system is nonlinear. So let's look at the nonlinear solution first.

The parameters:

a1 = {{-3, 2}, {-0.25, 1}};
a2 = {{-1.9, -0.4}, {-2.24, -4.7}};
b1 = {{0.25}, {1}};
b2 = {{-2.5}, {1}};
c = {{1, 0.5}, {0, 1}};
Subscript[ρ, 1] = (1 - Tanh[Subscript[x, 1][t]])/2;
Subscript[ρ, 2] = 1 - Subscript[ρ, 1];

Set up the nonlinear equations and obtain the solution.

xx = {Subscript[x, 1][t], Subscript[x, 2][t]};
eqns = Thread[{Subscript[x, 1]'[t], Subscript[x, 2]'[t]} == 
  Subscript[ρ, 1] (a1.xx + b1.{u[t]}) + Subscript[ρ, 2] (a2.xx +b2.{u[t]})];
ics = {Subscript[x, 1][0] == 0.1, Subscript[x, 2][0] == 0};
sols = NDSolve[
Join[eqns /. u[t] -> 10 Sin[t], ics], {Subscript[x, 1][t], 
Subscript[x, 2][t]}, {t, 0, 10}];
p = Plot[Evaluate[c.xx /. sols], {t, 0, 10}]

enter image description here

The solution of the linearized system.

StateSpaceModel[eqns, xx, u[t], c.xx, t];
OutputResponse[{%, {0.1, 0}}, 10 Sin[t], {t, 0, 10}];

Compare it with the nonlinear one.

Show[p, Plot[%, {t, 0, 10}, PlotStyle -> Dashed]]

enter image description here

share|improve this answer
    
Thanks a lot for your help. I am trying to understand all the code. Nevertheless the solution from Matlab and Mathematica are similar but quite different. Probably because in Matlab I am using Euler to solve the ODE. $\rho_i$ is a convex function, so in theory the LPV system is as the name says linear. It is possible to plot from this code the values of $\rho_i$ and also $x(t)$. Thanks for your help I am trying to change from matlab to mathematica. – user70012 Jun 14 '13 at 14:57
    
First, you might try ode23 or ode45 in Matlab -- Euler approximations often go wrong. Second, this is not a linear system since $\tanh(x_1)$ is not a rational function. "Convex" and "linear" are not at all the same thing, so you should expect to see different behaviors. – bill s Jun 14 '13 at 15:00
    
yes I understand, the idea of this kind of system is to have a set of LTI system [$A_ix(t)+B_iu(t)$] multiplied by a convex gain scheduling function (GSF) $\rho_i$, in this case a nonlinear function. In conclusion a set of LTIs multiplied by a nonlinear function. Where the GSF depend in the measurable (or unmeasurable) parameter varying $x(t)$. The advantage of this kind of system is that the overall representation is closer to the nonlinear behavior. Again thanks a lot for your help, I learn new things of Mathematica reading your code. – user70012 Jun 14 '13 at 15:15
    
For the nonlinear system, use Plot[Evaluate[vars /. sols], {t, 0, 10}] with vars = xx for the states and vars = {Subscript[[Rho], 1], Subscript[[Rho], 2]} for the parameters. For StateSpaceModel use StateResponse for the states, and change the outputs to obtain the parameters. Regarding nonlinear, I used it in the sense of not being linear (superposition and homogeneity fails.) – Suba Thomas Jun 14 '13 at 15:21
    
Your right. People in control also call this kind of systems as Takagi-sugeno. Thanks for your help. – user70012 Jun 14 '13 at 15:27

Update

As of Version 10 nonlinear control systems can be addressed by NonlinearStateSpaceModel. With this the answer given by @Suba Thomas can be obtained by:

a1 = {{-3, 2}, {-0.25, 1}};
a2 = {{-1.9, -0.4}, {-2.24, -4.7}};
b1 = {{0.25}, {1}};
b2 = {{-2.5}, {1}};
 c = {{1, 0.5}, {0, 1}};
ρ[1] = (1 - Tanh[x[1][t]])/2;
ρ[2] = 1 - ρ[1];

xx = { x[1][t], x[2][t] };
eqns = Thread[{x[1]'[t], x[2]'[t]} == ρ[1] (a1.xx + b1.{u[t]})
       + ρ[2] (a2.xx + b2.{u[t]})];
ics = {0.1, 0}; (* initial conditions for x[1][0], x[2][0] *)

(* Nonlinear State Space Model *)
nssm = NonlinearStateSpaceModel[ eqns, xx, u[t], c.xx, t ];
nssmResp = OutputResponse[ { nssm, ics }}, 10 Sin[t], {t, 0, 10} ];
nssmStyle = Sequence[ Directive[Blue], Directive[Red] ];

(* Linearized State Space Model *)
ssm = StateSpaceModel[eqns, xx, u[t], c.xx, t];
ssmResp = OutputResponse[ { ssm, ics }, 10 Sin[t], {t, 0, 10} ];
ssmStyle = Sequence[Directive[ Dashed, Blue], Directive[Dashed, Red]];

(* Plotting the Responses *)
legend = LineLegend[
    { Blue, Red, Directive[Dashed, Blue], Directive[Dashed, Red] },
    {"nonlinear g(x1)", "nonlinear g(x2)", "linearized g(x1)", "linearized g(x2)"}
];
Plot[ {ssmResp, nssmResp}, {t, 0, 10},
    PlotStyle -> {ssmStyle, nssmStyle}, 
    PlotLegends -> legend
]

OutputResponses

share|improve this answer
    
Note that I shun Subscript as the devil will holy water (repeating myself)... – gwr May 25 at 15:14
    
Yes, this is now the more natural way to obtain the solution after NonlinearStateSpaceModel came onto the scene. – Suba Thomas May 26 at 14:51
    
@SubaThomas I must admit that (not coming from engineering) I am still struggling with the question when to use something like ParametricNDSolve or the block based modeling of NonlinearStateSpaceModel. There seem to be ups and downs for both: WhenEvent, speed of simulation when varying a lot of parameters, state resetting for maybe Bayesian estimation ... Maybe that is a question worth asking in itself? – gwr May 26 at 15:03

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.