How do I make a plot with all three states [x1 x2 x3]' in SImulink?

I have made a model on people [x1 x2 x3]' influencing each other by influence matrix A.

How can I plot the changes in the states over k (discrete) and time (continuous)?

The initial value is x0 which has values for

x = [x1 x2 x3]'

with

x0 = [12 7.57 5]' 

And matrix A is a 3x3 matrix with unstable eigenvalues.

My discrete model looks like this:

And my continuous model ike this:

[Ameer Hamza]

Did you open the scope block? It should be showing three lines. As you know that the system is unstable, so they might not be easily visible due to the scaling issue.

[DB]

Thank you for the answer!

The scope reads the same as the dashboard scope with just one output plot, even when zoomed in.

It is also strange that the output is zero until the instability kicks in. This should not be the outcome.

[Ameer Hamza]

Can you attach your model file? Also, the signal is not zero before instability kicks in. It is just a matter of scale. During instability, there is an exponential growth in the signal. Therefore, the previous values are negligible as compared to new values. For example, inside the scope, you can use cursor measurements to check the values.

[DB]

I have attached my Matlab and simulink files for the discrete as the continuous time model. I hope you can help me!

[Paul]

Do the simulink models represent the actual equations you’re trying to solve? Making the initial state a constant input looks peculiar.

It sounds like the equations should be:

Xdot = Ac*X , X(0) = X0

or

X(k+1) = Ad*X(k), X(0) = X0

If so, you will need to modify your block diagrams.

[DB]

@Paul Yes that is what I intended! The input is constant as this is no ordinary control system. It is about consumption of people, and it starts at a constant value. I did not define xdot as I thought defining it in simulink would be sufficient. I will try it! I have Ac=Ad=A as it is constant. I will try x(0)=x0, I thought I got an error when putting in 0 as time or time step. Thank you for your answer!

[Paul]

So your continuous time equation is

Xdot = Ac*X + Bc*U, X(0) = 0, U = constant ?

Whatever your equation is, you can easiily represent it in Simuink with either a state space block with U as an input and X0 as the initial condition block parameter, or cobining an integrator and gains with X0 as the initial condition block parameter for the integrator.

Also, if the discrete time model is supposed to be an approximation to the continuous time model, then it's highly unliekly that Ac = Ad is correct.

[DB]

@Paul Thank you for your help.

The different thing of this model is that is has no control input, so I solely have an intitial state vector x0 and constant influence matrix A as inital point. I see now that B should be zeros(3) because of this. My discrete model now shows the desired behaviour, but it does not show the three states seperately but as the sum in one plot. How can I see them seperately?

The ultimate goal is to have a graph like the one shown below. The A should be constant as it applies to both continuous as discrete.

I have the continuous one corrected as seen below. But this does not show as promising results as the discrete one.

[Paul]

I still don't think the simulations are set up correctly based on the description of the problem. The initial state should not be modeled as a constant input. Note that in the continuous case, the initial states are all zero on the plot. The feedback loop around around the discrete time state space model is suspect as well, and I don't understand why a state space model is used for the discrete time case when only a delay is needed (analgous to the integrator in the continuous case).

Do you really need to do this in Simulink? The problem you've described can be easily implemented in Matlab alone. The Control System Toolbox can also be used.

Let's try to first get a handle on the problem you're trying to solve. Let's focus on the continuous case first. As I understand it, you want to solve the differential equation

xdot(t) = A*x(t),  x(0)= x0 = [12 7.57 5].'  

where A is a constant 3 x 3 matrix. To simplify this discussion, let's assume a diagonal A

A = diag([-1 -2 -3])

Because of the form of A, we know the solution to the differential equation

x1(t)=x0(1)*exp(-1*t)
x2(t)=x0(2)*exp(-2*t)
x3(t)=x0(3)*exp(-3*t)

We can compare this closed-form solution to our computed solution, which we develop as follows. In general, the solution to the given differential equation is

x(t)=e^(A*t)*x0

where e^(A*t) is the matrix exponential.

Let's put this all together and compare the computed solution to the close form solution.

x0 = [12 7.57 5].';
t = 0:.01:5;
A = diag([-1 -2 -3]);
x = nan(numel(t),3);
for ii = 1:numel(t), x(ii,:) = (expm(A*t(ii))*x0).'; end
xclosed = [x0(1)*exp(-1*t(:)) x0(2)*exp(-2*t(:)) x0(3)*exp(-3*t(:))];
plot(t,x,t(1:10:end),xclosed(1:10:end,:),'o'),grid

As must be the case, the solution starts from the initial conditions and decays to zero.

You can also get the same solution using the Control System Toolbox.

sys=ss(A,zeros(3,1),eye(3),zeros(3,1));  % Use eye(3) for the C-matrix so that each state is an output
y = initial(sys,x0,t);
plot(t,y,t(1:10:end),xclosed(1:10:end,:),'o'),grid

Either of these methods can be used with your particular A-matrix (which would be helpful to post if you can).

Is this the type of solution you're expecting? If it isn't, you'll need to provide more clarity on your problem. If it is, then the next step would be to use a similar approach for the discrete time case, and we can look at Simulink implementations if that's really needed.

[DB]

@Paul thank you for your response, I really appreciate it.

Simulink is not necessary, but was recommended to me. As long as I get the desired graph.

Unfortunately the first option does not give the desired result. Why did you include a matrix exponential? This function should be on normal time-scale, not logarithmic. The solution of the three states should converge, not nullify, as shown in the desired example graph.

The example A matrix is:

a11=0.95;   %x1 is stubborn because of high hedonic and gain goals
a12=0.0;    %x1 has no contact with x2
a13=0.05;   %x1 has a talk with x3 on energy usage
a21=0.05;   %x2 hears rumors of the high energy usage of x1
a22=0.7;    %x2 is easily influenced: not stubborn.
a23=0.25;   %x2 talks to x3
a31=0.02;   %x3 consumes some more because x1 does it too
a32=0.01;   %x3 consumes some more because x2 does it too
a33=0.97;   %x3 has high biospheric values

The secónd method does not show convergence either.

[DB]

I have a code that almost shows the desired continuous behavior. I will also need the discrete of it.

clear all
close all
a11=0.95;   %x1 is stubborn because of high hedonic and gain goals
a12=0.0;    %x1 has no contact with x2
a13=0.05;   %x1 has a talk with x3 on energy usage
a21=0.05;   %x2 hears rumors of the high energy usage of x1
a22=0.7;    %x2 is easily influenced: not stubborn.
a23=0.25;   %x2 talks to x3
a31=0.02;   %x3 consumes some more because x1 does it too
a32=0.01;   %x3 consumes some more because x2 does it too
a33=0.97;   %x3 has high biospheric values
A=[a11 a12 a13; a21 a22 a23; a31 a32 a33];
n=100;
x0 = [12 7.57 5].';
t = 0:0.01:5;
x(:,1)=A*x0;
for i=2:1:n
    x(:,i) = A*x(:,i-1);
end
t=1:1:n;
figure;
plot(t,x);

[Paul]

I must not understand what you mean by "discrete" and "continuous." You said that this code shows the "continuous" behavior, but the solution you posted here is the solution to the dicrete time difference equation

x(k+1)=A*x(k); x(0)= x0

which is exactly the equation that is implemented in the code

x(:,1)=A*x0;
for i=2:1:n
    x(:,i) = A*x(:,i-1);
end

Whether or not the states converge to nonzero values, diverge, converge to zero, etc. depends on the particular A matrix and the initial condition of the state. In this particular case the eigenvalues of the A matrix are

>> eig(A)

ans =

   6.9126e-01
   9.2874e-01
   1.0000e+00

Because the first two are inside the unit circle and the last one is on the unit circle, the solution will converge to a non-zero value, as shown above. If you want to see a different behavior, try

a11 = 0.8;

Even if this value may not be physically realistic for your use case, it will show you a different behavior.

[DB]

The graph shows continuous behaviour as it does not have the discrete "stairs-like" plot as is shown in the example graph. It look more like the continuous part of the example graph. The three states do converge to the same value as desired as values of A are all less than 1.

[Paul]

Do you just want to make a plot that has that staircase look?

stairs(t,x.')

[DB]

That plot function always gives an error message at my version. But the stairs plot is another way of plotting. I want to have a discrete prediction model with continuous approximation.

[Paul]

I'm not sure how much more (if any) help I can offer. I'm used to doing the opposite, i.e., starting with a continuous time model (differential equation) and then coming up with a discrete time approximation (difference equation). You seem to want to go the other way. One could come up with a differential equation that is approximated by the difference equation, but the outputs of those equations in response to the initial condition will not look like the plots you've posted.

The plot you posted shows a group of three (?) red lines and two blue lines. As best I can understand based on the problem description, the blue lines cannot be initial condition responses from discrete and continuous models where one model is the approximation of the other. One can come up with a model with two outputs that have the same appearance as the blue curves, but I wouldn't call that model a discrete model with a continuous approximation.

Sorry I couldn't be of more help.

[DB]

@Paul Thank you very much for your help. I do acknowledge I have a vague problem and I'll work on it further. I think your help has made me understand it a whole lot better.