State Space with Disturbance

I have this system:

Xdot = AX +Bu +Fd,

where A,B and F are the system, input distribution and disturbance distribution matrices respectively and X, U and d are the state, control and disturbance vectors respectively.

i have all the matrices data. i want to know how to use matlab to plot the state variables for a step disturbance.

where i am finding problem is because the standard form of state space in matlab is :Xdot = AX +Bu without the distubance vector.

so how do i input my system in Matlab???

Answer 1

[Paul]

Combine d and u into a single input vector:

Xdot = A*X + [B F] * [u;d]

y = C*X + [D 0]*[u;d] % assuming the disurbance doesn't feed directly into output

sys = ss(A,[B F],C,[D zeros(size(C,1),size(F,2))])

[Manuel Valerio Antoni]

This answer is very interesting for me too. The only error I get with that, is this: Dimensions of arrays being concatenated are not consistent. In fact I have a 12x4 B matrix and a 12x6 F matrix, with of course 4 inputs and 6 disturbances. Is there a simple way to solve this?

[Paul]

Can you post an example? Seems like it should work:

A = rand(12);
B = rand(12,4);
F = rand(12,6);
C = rand(1,12);
D = rand(1,4);
plant = ss(A,B,C,D);
augmentedplant = ss(A,[B F],C,[D zeros(size(C,1),size(F,2))]);
disp('ran without error')

ran without error

[Manuel Valerio Antoni]

Pardon, my error. I didnt' specify a proper dimension for D (just set it =0; previously). However, how can I implement this new definition in Simulink if the "State-space" block is defined in the canonical form and it's not possible to specify the disturbances matrix?

[Paul]

Not sure what you mean by canonical form.

Assuming you have A, B, C, D, and F defined in a workspace ....

In the State Space block enter the dialog parameters as:

A -> A

B -> [B F]

C -> C

D -> [D zeros(size(C,1),size(F,2))]

The input to the State Space block should be formed by combining the signals u and d with a mux block (u on top, d on the bottom).

[Mohammed Yakoob]

Hello Paul and thanks for this explanation but my question is how can I deal with this disturbance if using the LQR controller?

[John]

@Paul Just to check, does your augmenting vectors (where [u;d] are the inputs) allow independent control of u and d using lsim()? I ask, because I tried the below, and am seeing an unexpected (to me) result where an input disturbance of 1 does not give a position shift of 1.

For example, assume some basic 2nd order system for the open-loop SS model (m = 1e-5, k = 62534, b = 0.0016)

openLoopSsSys = ss(A, [B F], C, D)

with A = [0 1; -6.3e4 -160], B = [0; 6.3e4]

x1 = position

x2 = velocity

Assume the lsim() time is 0 to 1 seconds, and the input parameter to lsim() is [u(t); d(t)] , instead of just the standard u(t), as you proposed:

lsim(openLoopSsSys, [u; d]).

Matlab seems to correctly simulate a closed loop response to a reference input, and closed-loop response to a disturbance, if the input vector is defined as (for example) below

u = 1 in t = [0, 1] (Ie all the simulation time)

d = -1 if t > 0.75, d = 0 otherwise;

I do see some change, but not as expected.

When I added d = -1 as above, there was virtually no x1 change output. The "d" input must be much different -1 (-400), even though the mapping above should be:

x1_dot (position change) = x2 + d + u

I thought an input d = -1 would shift x1 (position) by -1, instead of needing to be (in this example) -400.

Why is that? Is lsim() somehow interpreting this differently than i thought, or is my understanding incomplete on how d maps to x1?

Thanks.

[Paul]

Yes, the inputs to lsim can be whatever you want. Of course, the system input to lsim has to be defined properly for the system you want to represent, but that's a different problem. For example, suppose we have a plant model

xdot = A * x + B * u + F * d

y = C * x + D * u % assuming d does not directly feed y to simplify things a bit

Using state feedback:

u = r - K * x

xdot = A * x + B * (r - K * x) + F * d = (A - B*K) * x + B * r + F * d

y = C * x + D * (r - K * x) = (C - D*K) * x + D * r

sys = ss(A - B*K, [B F], C - D*K, [D zeros(size(C,1),size(F,2))];
t = 0:.01:1;
u = 1 + 0*t;
d = -1.*(0.5 <= t & t <= 0.75);
y = lsim(sys,[u(:) d(:)]);

[Paul]

Hi John,

Looks like we were both typing at the same time. Trying to reconstruct the example

m = 1e-5; k = 62534; b = 0.0016;
h = minreal(tf(1,[m b k]))  % tf from force to position

h =
 
           1e05
  ----------------------
  s^2 + 160 s + 6.253e09
 
Continuous-time transfer function.

sys = ss([0 1;-k/m -b/m],[0;1/m],[1 0],0)

sys =
 
  A = 
               x1          x2
   x1           0           1
   x2  -6.253e+09        -160
 
  B = 
          u1
   x1      0
   x2  1e+05
 
  C = 
       x1  x2
   y1   1   0
 
  D = 
       u1
   y1   0
 
Continuous-time state-space model.

tf(sys)  % verification

ans =
 
           1e05
  ----------------------
  s^2 + 160 s + 6.253e09
 
Continuous-time transfer function.

That doesn't match your system, the B-matrix is different. Maybe your model is different than what I assumed?

Also, based on the model posted

format short e
A = [0 1; -6.3e9 -160], B = [0; 6.3e4]

A = 2×2

1.0e+00 * 0 1.0000e+00 -6.3000e+09 -1.6000e+02

B = 2×1

0 63000

The equation for x1dot is (after adding the disturbance into the problem): x1dot = x2 + 6.3e4 * u - d

So u will have a much, much bigger effect on x1dot than d.

If you'd like to continue the discussion, I suggest posting a complete code, show the results from running that code, and describe what about those results are of concern.

[John]

Thanks @Paul. You're right, that was confusing. Let me start over and be more specific with another related question, that's more fundamental to my questions. As you requested a test file is attached.

Here's an open-loop 2nd order system with a disturbance vector added per your recommendation. Since it's a basic 2nd-order spring-mass-damper system, and the state space was set up in the standard way, the output states are x1=position, x2=velocity. I'm plotting the lsim() output states x in response to a constant reference step at t=0 (the sys is normalized), and dist step at t = 0.25.

When plotting all state variables, the position changes as expected, and the velocity decays to 0 as expected after settling after a reference command.

But the velocity stays constant after a disturbance (force) input. Presumable, since the position stabilizes (forces in equilibrium), velocity should be 0.

As Matlab treating this input vector correctly?

(The code should generate this plot as well)

[Paul]

Hi John,

The short answer is that I think the intent was to add in the disturbanc as a force, i.e., at the acceleration level, but it was actually added in at the velocity level.

>> sysOlWithDistDc

sysOlWithDistDc =
 
  A = 
                  x      x dot
   x              0          1
   x dot  -6.25e+04       -150
 
  B = 
            r(ref)   F(dist)
   x             0     416.7
   x dot  6.25e+04         0
 
  C = 
            x  x dot
   pos      1      0
 
  D = 
         r(ref)  F(dist)
   pos        0        0
 
Continuous-time state-space model.

The equation for x1dot is: x1dot = x2 + 416.666666 * Fdist. So, with Fdist = -1 the steady state value of x2 will be 416.6666666, which is exactly what is seen. More generally, for an input of r = 1 and d = -1, we know this system will have an equilibrium point with with states

>> -sysOlWithDistDc.A\(sysOlWithDistDc.B*[1;-1])

ans =

                         0
     4.166666666666666e+02

So the equilibrium position is x1 = 0, also as seen.

If the disturbance was intended to be a force as indicated in the comments, then change this line

%% Add open-loop disturbance
% Disturbance input is a force input (not position)
% F = [1; 0];
F = [0;1]; % force sums into x2dot !

and you'll get a reponse you're expecting with x1 = constant and x2dot = 0 in the steady state.

Also, the code is using a variable called Ctf, which isn't defined, but that is easily fixed.

I'm not quite sure what the intent is of all the dcgain scaling .... but that doesn't matter for this analysis, other than how it affects the steady state values.

[John]

Thank you, @Paul. That starts to clarify. However, if this is the case, how can x remain constant, if xdot is non-zero constant due to a velcity disturbance input? The incompatability of x and xdot is what's most confusing.

Regarding the scaling: it's just to get a unit output for a unit step input. I'm just normalizing for a nice plot that is easy to digest, akin to a gain in a physical system calibration of command input to motion.

[Paul]

Because in the model x2 is not the derivative of x1. If it were, the model would have the equation

x1dot = x2

But that's not the equation in the model. The equation in the model for x1dot is

x1dot = x2 + 416.66666 * Fdist

Therefore, in the steady state it must be true that x1dot = 0, which in turn implies that in steady state we must have x2 = -416.6666 * Fdist, which is exactly what happens. Because x2 is not the derivative of x1.

[John]

1) Thanks @Paul, I see. But then since x2 = k*F - x1dot, how is Fdist a velocity disturbance input to the original spring/mass/damper system?

It's just connecting a disturbancec only to specific states, and breaks the original SS equation set that describes the physical system.

I'm looking for disturbance inputs that stimulate physical systems that behave like one.

Ie with a velocity input, position should also change (x1=ramp, x2=Const).

With a force input, equilibrium should be the result (x1=Const, x2=0).

So isn't this a non-physical system with the dist spec-d as above? How would I spec a vel dist input that affects the spring/mass/damper?

Or what am I missing? :)

2) How is the Matlab native K disturbance matrix used in sys() relate to this strategy of augmenting a disturbance to SS? K seems a mapping, presumably to take e(t) to the appropriate state variable (or combinations).

"sys = idss(A,B,C,D,K) specifies a disturbance matrix K.

"A state-space model of a system with input vector u, output vector y, and disturbance e takes the following form in continuous time:

dx(t)dty(t) = Ax(t)+Bu(t)+Ke(t)"

(a) how does K support multiple inputs (eg ref and dist)

(b) how is K used concretely? The documentation has no example for an input e(t).

Is K even suggested for use? :D

[Paul]

1). The equation is x2 = x1dot - K*F. Because the system is stable, we know that x1dot will be zero in response to constant inputs. So, in steady state we have x2 = -K*F, which is exactly what happens.

At this point, I'm not sure what you're trying to do. I only know what you've done. The model you've defined is

x1dot = x2 + 416.7*Fdist

x2dot = -6.25e34*x1 - 150*x2 + 6.25e4*rref

and that's the model that lsim simulates.

If those equations don't model the system as intended, then the equations need to change. If, as I suspect, the intent is that x1 is position, x2 is velocity, and Fdist is a force, then the equations need to be changed such that Fdist is summed with the rest of the forces in the system.

x1dot = x2

x2dot = 1/m*(-k*x1 - b*x2 + rref + Fdist)

assuming that rref is force (and before any of that dcgain scaling). These equations represent a physical system with two externally applied force inputs, rref and Fdist.

2. I'm only vaguely familiar with the System ID Toolbox and so can't comment on usage of idss with non-zero k.

[John]

Thanks @Paul; I appreciate your help. I'll try to be more clear; apologies for being confusing:)

"At this point, I'm not sure what you're trying to do."

I'll try to restate the question more clearly, and ask a concrete question at the end.

I'm modeling a basic open-loop 2nd order spring-mass-damper system with a disturbance vector input, to model disturbances (F augmentation per your recommendation). The state space was set up in the standard way as in my post above; so the output states are x1=position, x2=velocity=x1dot. The physical motion as a result of the disturbance inputs, must always be representative of how this physical system would actually behave.

Using ss() and lsim(), I'd like to both model

a) a force F to create a position disturbance, eg a step F, and

b) a velocity disturbance -- needs a changing F, I think -- eg a nonzero vel in steady state (after system has been accelerated). It could be an infinitely accelerating frame so the dF is constant and leads to constant velocity (and thus ramping position change). Or it could be eg a shaking table the system's reference frame is attached to, so though the disturbance vel mean could be zero still, but steady state is still not zero vel.

I'd like to plot lsim() output states x in response to (a) and (b). Let's start with velocity since that's what the .m file has, and since you mentioned that I already had a "velocity disturbance". (Can deal with F another time). A true velocity disturbance should cause the system's position (x1) to change continuously.

In a physical system, this nonzero vel (x2=x1dot) would only be caused be a changing F I think; a steady-state F would be a position disturbance (or an intentional reference input), and cause equilibrium and x1dot-->0. Per above, this nonzero constant vel could be an accelerating system, for example, or if non-constant vel, it could be a shaking system.

Questions: For

[B, F]; F = [1; 0];

What do you mean by "velocity disturbance" augmentation in the context of this physical system? Do you just mean a disturbance into the measurement? If so, why is it inputted to the state vector?

It doesn't seem like it's what I'm hoping for ("someone disturbs the spring/mass/damper system to cause it to move with some velocity for a long time, ie not a position step change"), since velocity disturbance steps are impossible for a system with mass; only the requested dF is possible.

Best,

John

[Paul]

Perhaps the problem is that the system that is represented by the state space model is not the system that you actually want to represent. If you want x1 to be the position of the mass relative to a fixed reference frame, and if you want x2 to be the velocity of the mass relative to the same reference frame, the model must have the equation

x1dot = x2

Since that's not the equation in your model, either the model is incorrect, or you want to model a different system. As best I can tell, that's the equation you want in your model.

The second equation should just be x2dot = (sum of forces)/mass, or

x2dot = (-k*x1 - b*x2 + u + d)/m,

where u (what you called rref) is a force input, perhaps one that will ultimately be determined from feedback control, and d is force disturbance that is some function of time applied by Mother Nature. So the second equation is

x2dot = -k/m *x1 - b/m * x2 + u/m + d/m

We can write these equations as

xdot = [0 1; -k/m -b/m] * x + [0;1/m] * u + [0;1/m] * d = A * x + B * u + F *d = A * x + [B F] * [u;d]

If the ouputs are just x1 and x2, then we have

y = I * x = C *x.

Putting in the numbers:

zeta = 0.3;
wn = 250;
m = 1e-5;
b = 2*zeta*wn*m;
k = wn^2*m;

A = [0 1;-k/m -b/m];
B = [0;1/m];
F = [0;1/m];
C = eye(2);
D = zeros(2,2);
sys = ss(A,[B F],C,D,'StateName',{'pos' 'vel'},'InputName',{'u' 'd'},'OutputName',{'pos' 'vel'})

sys =
 
  A = 
              pos        vel
   pos          0          1
   vel  -6.25e+04       -150
 
  B = 
            u      d
   pos      0      0
   vel  1e+05  1e+05
 
  C = 
        pos  vel
   pos    1    0
   vel    0    1
 
  D = 
        u  d
   pos  0  0
   vel  0  0
 
Continuous-time state-space model.

Now simulate the response of the system to constant input u = 1 for all time, and constant disturbance d = 5 starting at t = 0.5

t = (0:1e-3:1).';
u = 0*t + 1;
d = 5*(t > 0.5);
y = lsim(sys,[u d],t);
figure
plot(t,y(:,1)),ylabel('pos')

plot(t,y(:,2)),ylabel('vel')

As expected, we see a lightly damped response, the steady state position is constant, and the steady state velocity is zero.

Now change the disturbance input so that it's a ramp force instead of a step

d = 50*(t-0.5).*(t > 0.5);
y = lsim(sys,[u d],t);
figure;
plot(t,y(:,1)),ylabel('pos')

The position goes to a constant until the ramp force disturbance hits, and then the position ramps up as well.

Plot the velocity and zoom in.

plot(t,y(:,2)),ylabel('vel')
xlim([0.4 1])

As expected, steady state velocity is zero for t < 0.5, and then goes to a constant after the disturbance hits, which should be the slope of the position (after the transient due to the disturbance dies out).

These results are all as expected based on the model and are consistent with expectations based on the physics of the problem.

Now, if we change the state-space realization of sys such that F = [1;0], we are no longer representing the spring/damper system we started with. It's now the state-space model of a different system. That system might not even represent a physical system, but math and Matlab don't care, even if you don't change the StateName or InputName. The math is allowable and Matlab will compute the response for that ss model, whatever that response may be.

BTW, I don't think I ever used the term "veclocity disturbance," at least I didn't mean to. I did say " it was actually added in at the velocity level." Here, "velocity level" just means that the disturbance was added into the x1dot equation, as opposed to adding it to the x2dot equation.

[John]

Thanks @Paul.

1) "Perhaps the problem is that the system that is represented by the state space model is not the system that you actually want to represent"

Huh...could you clarify?

I thought the system in the .m file i shared, with x1dot = x2, is the same as your post:

A = [0 1; -k/m -b/m];

B = [0; 1/m];

is from the .m file, and

A = [0 1;-k/m -b/m];

B = [0;1];

is from your code above (I think you might have forgotten to divide by m in B, but that doesn't impact anything beyond scaling).

What's the difference? :D Sorry, i'm guessing it's obvious, but I'm missing it...

2) "BTW, I don't think I ever used the term "veclocity disturbance,"

Ah, quite so -- you're absolutely right. Apologies for misrepresenting!

3) "Now, if we change the state-space realization of sys such that F = [1;0], we are no longer representing the spring/damper system we started with." "The math is allowable and Matlab will compute the response for that ss model, whatever that response may be."

Got it -- that makes sense. I'd overlooked that the SS formulation specifically only makes physical sense for a force input at the acceleration level (F=ma), ie x2dot=x1: level.

[Paul]

1). I was referring to the state space model called sysOlWithDistDc, which included the disturbance as an adiditonal input that was added into the x1dot equation, whereas here the disturbance is added in to the x2dot equation. I did forget to divide by m in B and F. I went back and fixed that.

[John]

Thanks @Paul, that all makes sense. Thanks for the explanations and clarifications; this was quite informative.

[Wei Zhao]

Thank both of you! Very happy to find this solution.

[jana nassereddine]

but in this way the disturbances is considered as a manipulated variable, how can we add the disturbance as an input to model predictive control, and with matrix related to the state of the plant model?

[Paul]

Hi jana,

The state space model can have any number of inputs as required to represent the dyanmics of the system under consideration. Some of those inputs will be control inputs, which is what I think you mean by "manipuated variable," that we can influence with compensation, and other inputs will be disturbances that we have no control over. In the Control System Toolbox representations, both types of inputs are lumped together in a single input vector in the model.

Can you show a simple example for what you're concerned about with model predictive control?

[Tomo]

Hello @Paul,

If you model it this way like you said,

Xdot = A*X + [B F] * [u;d]

Does that mean to determine reachability we would compute W = ( B_{new}, AB_{new},...) where B_{new}=[B F].

Thanks.

[Paul]

If reachability is to be assessed based only on the control input to the system (as is typical), then the reachability matrix would be computed from A and B.