You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
State Space with Disturbance
254 views (last 30 days)
Show older comments
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???
Answers (1)
Paul
on 21 Jan 2021
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))])
23 Comments
Manuel Valerio Antoni
on 27 Sep 2021
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
on 27 Sep 2021
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
on 28 Sep 2021
Edited: Manuel Valerio Antoni
on 28 Sep 2021
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
on 28 Sep 2021
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
on 4 Apr 2022
Hello Paul and thanks for this explanation but my question is how can I deal with this disturbance if using the LQR controller?
Sam Chak
on 4 Apr 2022
Are you asking if the LQR can reject the disturbance, or, how to reject the disturbance using MATLAB Control System Toolbox if the LQR is already used?
I'd suggest you to post a new question with the dynamical system and disturbance. Then, tag @Paul like this to request for his assistance or explanation in the disturbance rejection-based control design.
John
on 24 Feb 2023
Edited: John
on 24 Feb 2023
@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
on 24 Feb 2023
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
on 24 Feb 2023
Edited: Paul
on 24 Feb 2023
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
on 25 Feb 2023
Edited: John
on 25 Feb 2023
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
on 25 Feb 2023
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
on 26 Feb 2023
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
on 26 Feb 2023
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
on 26 Feb 2023
Edited: John
on 26 Feb 2023
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).
"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
on 27 Feb 2023
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
on 27 Feb 2023
Edited: John
on 27 Feb 2023
"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
on 27 Feb 2023
Edited: Paul
on 27 Feb 2023
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
on 27 Feb 2023
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
on 27 Feb 2023
Edited: Paul
on 27 Feb 2023
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.
jana nassereddine
on 6 Apr 2023
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
on 6 Apr 2023
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?
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)