How to add disturbance feedforward to LQI()

I'm seeing an issue where adding disturbance FF to an LQI system is not causing any performance increase -- even though in my Simulink results, disturbance FF by itself yields a much faster response than FB. It looks like LQI (by its nature of state feedback that doesn't generate a traiditional error signal) is treating the dist FF input as a disturbance, and canceling it out.

1) How can I use LQI() to add disturbance feedforward (FF), if i had a sensor that measures the disturbance w acting on sys?

Image from Matlab LQI() link:

2) Normally, it's easy to add disturbance FF to classical controllers like PID, since the FB and FF loops are decoupled; so the process is to just make a new u: u = u_FB + u_FF, and adding FF will just increase the performance and generally not impact FB stability.

Since LQI uses state feedback (not error-based feedback) for all parameters except y, any u_FF injected to sys will cause LQI (the FB loop) to try and cancel out the FF contribution, since u_FF looks like an unwanted disturbance that causes LQI to not track according to the target trajectory set by the designer.

So, adding u = u_FB+u_FF causes the exact same response as u = u_FB, ie no performance gain.

3) I have several Simulink sims of this i can share if needed, but I'll start by asking the basic conceptual question first :)

My sims show FF-only acting about 4x as quickly as FB-only. And FF+FB acts exactly as quickly as FB only, since u_FB also then adds a mirrored component (-1*u_FF) to cancel it out, since FB is trying to maintain its designed trajectory.

FF only: settling time T given a dist. u = u_FF

FB only: settling time 4*T given a dist. u = u_FB

FB+FFonly: settling time 4*T given a dist (same as FB only). u = u_FB' + u_FF, where u_FB' = u_FB - u_FF. So u = u_FB.

4) Would something like the LQG Track method support feedforward, since it explicitly includes w? If so, how would I use LQG Track?

https://www.mathworks.com/help/control/ref/ss.lqgtrack.html

Answer 1

[Paul]

Hi John

I don't udnerstand your observations absent additional information. If the control input includes the negative of the disturbance (assumed to be perfectly measured), then it seems to me that the disturbance should have no impact on the output.

Here's a simple plant model, with a control input (u), a disturbance input (w), an output measurement from the plant (y), and an output that's a measurement of the disturbance (wmeas).

plant = ss(-1,[1 1],[1 ;0],[0 0;0 1],'InputName',{'u' 'w'},'OutputName',{'y' 'wmeas'})

plant =
 
  A = 
       x1
   x1  -1
 
  B = 
       u  w
   x1  1  1
 
  C = 
          x1
   y       1
   wmeas   0
 
  D = 
          u  w
   y      0  0
   wmeas  0  1
 
Continuous-time state-space model.

State feedback gain to place the closed loop pole at s = -4.

kfb = 3;

Define the controller as u = kfb*(r - y) - B(2)*wmeas, i.e., augment the feedback control with the negative of the disturbance input to the plant.

controller = ss([kfb -plant.b(2)],'InputName',{'e' 'wmeas'},'OutputName','u')

controller =
 
  D = 
          e  wmeas
   u      3     -1
 
Static gain.

Form closed loop system

sysc = connect(plant,controller,sumblk('e=r-y'),{'r' 'w'},{'y'})

sysc =
 
  A = 
       x1
   x1  -4
 
  B = 
       r  w
   x1  3  0
 
  C = 
      x1
   y   1
 
  D = 
      r  w
   y  0  0
 
Continuous-time state-space model.

We see by inspection that w has no impact on the system because it's multiplied by zero in the B and D matrices.

[John]

Thanks @Paul

1) "If the control input includes the negative of the disturbance (assumed to be perfectly measured), then it seems to me that the disturbance should have no impact on the output."

Maybe i wasn't clear enough and gave the wrong impression :) so here are more details. The control input is not the negative of the disturbance; It's that the FB cmd includes the negative of the command signal from FF (u_FF), so FB is canceling out FF, in order to maintain it's own disturbance-rejection profile:

u = u_FB' + u_FF, where u_FB' = u_FB - u_FF (this is not desired at all -- and is just showing that when FF is active, u_FB cancels out f_FF). So u = u_FB only, which is not expected or desired.

Say there's a disturbance at the output y, with dist FF in the manner shown below.

This is a stock image so ignore the rest of the loop eg Gc; i'm only showing this to show where dist FF enters the plant via a command. In my case Gc is LQI so the FB isnot all from E, since some are from direct state feedback.

2) In the plots below, i show the dist rejection on the top plot with FF only, and also FF + FB. As can be seen, FF can be very fast, but FB+FF is way slower. In fact, the FB + FF response is exactly the same as the FB-only, so i'm not plotting FB-only.

The bottom plot shows the FF+FB command int he FF+FB setup, where u_FB shows a -1*u_FF contribution, which cancel out FF (which i don't want, but that's what occurs). The cursor isn't showing anything interesting, i just forgot to turn it off.

3) Here's the model (set to FB+FF via the switches), and isn't made to be presentable, but hopefully it's clear...

I'm not showing what the Gff function is (it's basically inverse plant dynamics), and it's fed in from the top at u_FF. The dist_input at the output y is a step input.

Note: I have two such models running in parallel to make the plot above that shows both FF-only and and FB+FF, but i'm only showing one of the models; the other is identical but with differently-set switches.

As can be seen, the FB controller sees u_FF as an unwanted disturbance, and "correctly" rejects all the disturbances (as it sees it...though only dist_input is the true disturbance) to get a response that would occur from FB-only.

By the way, I'm using the paramterized method for LQI, to set the FB response to be a specific settling time, zeta, etc., and this FB+FF (and FB-only) response is perfectly matching that. So FB is working hard to reject the FF cmd.

4) I'm trying to understand how to tell LQI that u_FF is not an actual disturbance, but is helping. Ie FB shouldn't cancel it out, and should allow FF's contribution, which will help settle faster than the LQI settling time spec would be -- and this would be a good thing since dist rejection would be higher.

In other words: FB should allow for faster settling than the parameterized LQI spec, when FF is active.

The ideal is: FF works much faster, so allow FF to do it's thing, so that when FB contributes, it has less of a disturbance to act on with its predefined dynamics.

Maybe LQI Track does that? But I don't believe it does when i tried.

5) I'm also confused at how FB can reject FF perfectly, since that seems to indicate that FB would also be able to reject with FF's abilities. But that doesn't seem to be the case, ie i cannot set FB to reject the dist in full with FF's performance (the blue curve), and yet somehow FB can fully reject FF's contributions.

That seems to imply that there's a difference between canceling out dynamics that are just starting (ie from FF), vs being able to fully drive to a complete dist rejection (ie FB was the only contributor trying to match the blue curve). But i wouldn't expect the above being consistent with an LTI system (which the plant is: it's a 2nd-O LTI system)

[Paul]

Let me make sure I understand correctly.

The feedback signal to the integral control (ignoring the transport delay and noise inputs) is (in the frequency domain):

Ymeas(s) = Y(s) + D(s)

The control input to the plant is (in the frequency domain):

U(s) = -K*X(s) - Ki*E(s)/s + Gff(s)*D(s)

And the observation is that, for a step disturbance d(t) = step(t), the output y(t) (which is not the same as ymeas(t)) with Gff(s) = 0 (i.e., the FB-only case) is the exact same as the output y(t) with Gff(s) being whatever you've defined (i.e., the FF+FB case)?

[John]

Yes, that's correct.

And when Gff(s) !=0 and LQI is disabled (ie only FF is active), the response is much faster than FB alone. I'm unable to add both FF and FB together to gain higher performance, since FB is actively canceling the FF command when FF is enabled (see above plot on cmd u, for how FB includes a -1*u_FF), since it seems FB wants its own specific time constant for responses; see below for more on this.

PID+FF works as expected for FB+FF, I think since PID all works on feedback, so lower error arising from assistance by FF, is strictly better as far as the PID system is concerned. But with LQI's state feedback whose commands are not all mapped through an error feedback, it seems that LQI is actively regulating its response to be what i've specified, ie

H_lqi = [2 * zeta * wnTarget_rPs, 1, -wnTarget_rPs^2];
Q_lqi = H_lqi.'* H_lqi;
R_lqi = 1e-3;
[K_lqi, S1, P1] = lqi(plant, Q_lqi, R_lqi)

My theory is that LQI is trying to maintain wnTarget_rPs with zeta damping in response to any disturbance, and so cancels out u_FF, since u_FF causes a response that's much faster than the LQI spec -- aka FF looks like a disturbance as well, so FB tries to remove it. I can't see a way to let LQI know that FF is not a disturbance, but is a good input that should be left alone.

However, if i understand right what LQI is trying to do (and I might not be), FF causing a faster disturbance correction would fundamentally undermine an LQI who's trying to maintain specific response dynamics. Ie, somehow LQI needs to be okay with a response that's better than what it's been speced to achieve, when FF is active.

I tried feeding the FF signale through a plant model, and then into the setpoint reference, so that LQI sees the FF input as an appropriate trajectory. But (among other issues, including plant drift of that separate model path), LQI still regulates setpoint tracking through the integrator, so FF is again slowed to the FB dynamics.

[Paul]

Can you upload Klqi2, Plant2, and Gff in .mat file? Use the paperclip icon in the Insert menu.

[John]

Sure, attached. I also included a test slx.

In this example, arbitrary control efforts are possible, but that's okay even though it's nonphysical. I'm more trying to understand how to get a Simulink model to have FF and FB work together with LQI().

[Paul]

Load the data

load FF_check

Augment the plant to include the states as outputs and rename the output variables for easier typing (though perhaps less clarity)

Note that y = x1, which we'll come back to below.

plant2.C = [plant2.C;eye(2)];
plant2.OutputName = {'y' 'x1' 'x2'};

Feedforward control

gff2.InputName = 'd';
gff2.OutputName = 'uffin';

Add disturbance for measured output for feedback to the integrator. As done in the Simulink diagram, the disturbance is added to y, but is NOT added to x1, even though y and x1 are the same physical quantity. Make sure that's what you want.

sensor = sumblk('ymeas = y + d');

Feedback gains

Ki     = ss(K_lqi2(3),'InputName','ei','OutputName','ui');
Kstate = ss(K_lqi2(1:2),'InputName',{'x1' 'x2'},'Outputname','ustate');
% use this to include disturbance on "state" feedback as well.
%Kstate = ss(K_lqi2(1:2),'InputName',{'ymeas' 'x2'},'Outputname','ustate');

Switching gains

Kfb    = tf(realp('Kfb',1),'InputName','ufbin','OutputName','ufb');
Kff    = tf(realp('Kff',1),'InputName','uffin','OutputName','uff');

Form closed-loop system

sys = connect(...
    sumblk('e = r - ymeas'), ...
    tf(1,[1 0],'InputName','e','OutputName','ei'), ...
    Ki, ...
    Kstate, ...
    plant2, ...
    sensor, ...
    gff2, ...
    sumblk('ufbin = -ui - ustate'), ...
    sumblk('u     = uff + ufb'), ...
    Kfb, ...
    Kff, ...
    {'r' 'd'}, ...
    {'y' 'ymeas'});

Verify that the step response from d to ymeas is essentially the same using feedback with or without gff2

figure
step( ...
    sampleBlock(sys('ymeas','d'),'Kfb',1,'Kff',1), ...
    sampleBlock(sys('ymeas','d'),'Kfb',1,'Kff',0));

I believe this result can be explained as follows

Let P(s) be the transfer function from u to ang_vel

P = plant2('x2','u');

Break out the feedback gains

Ki     = -K_lqi2(3);
Kvel   =  K_lqi2(2);
Ktheta = K_lqi2(1);

Rename feedforward term for simplicity

F  = gff2;F.OutputName = '';

The closed-loop transfer function from d to ymeas is then

s = tf('s');
H = (F*P/s + 1 + Kvel*P + Ktheta*P/s)/(1 + Ki*P/s^2 + Ktheta*P/s + Kvel*P);
H = minreal(H);

15 states removed.

Verify by including on the step response plot

figure
step( ...
    sampleBlock(sys('ymeas','d'),'Kfb',1,'Kff',1), ...
    sampleBlock(sys('ymeas','d'),'Kfb',1,'Kff',0), ...
    H);

It seems to me that the strategy for selecting F(s) was to ensure that F(s)*P(s)/s + 1 is approximately zero at low frequency. However, that only takes care of two terms in the numerator of H(s).

Bode plot of all four terms in the numerator

figure
bode(F*P/s,tf(1),Kvel*P,Ktheta*P/s,{10, 1e7})
legend('F*P/s','1','Kvel*P','Ktheta*P/s')

We see that the F(s)*P(s)/s and 1 effectively cancel each other to roughly 1000 rad/sec. However, the numerator of H(s) is dominated by Kvel*P(s) and Ktheta*P(s)/s. So using F(s) to cancel the 1 doesn't really do anything. If ymeas is used as the feedback signal to Ktheta (or in the commented out line of Kstate), then Ktheta*P(s)/s won't show up in the numerator of H(s), but the Kvel*P(s) term will still dominate.

So, in the current form of the system (d only affecting feedback to integrator), it seems like F(s) should be selected such that it is realizable and satisfies

F(s)*P(s)/s + 1 + Kvel*P(s) + Ktheta*P(s)/s = 0 at low frequency.

But I didn't try to define F(s) this way and test it.

[John]

Thanks @Paul, that's useful. I didn't know that signal paths could be turned on and off as you showed in sampleBlock. Very handy :)

"As done in the Simulink diagram, the disturbance is added to y, but is NOT added to x1, even though y and x1 are the same physical quantity. Make sure that's what you want."

Perhaps that's what i meant...but not sure now. An example: say the system is a flying platform, trying to maintain some alignment to a target (ie a regulator), with the target position being measured by y; x1 and x2 are the states. The disturbance is the target suddenly shifting away, so although the platform has the same attitude (x1, x2), y changes. The goal is to shift the platform to aim at the new target. This is not a setpoint change, since the change in target shows up as a change in the error, so the flying system is still acting as a regulator.

So, i didn't assume they were the same physical quantity, but if they are, isn't what I have the proper simulink representation of a system who's output measurement is being disturbed or having noise added to it? Eg this LQGTrack() link. That should only affect y, and not x1.

"So using F(s) to cancel the 1 doesn't really do anything. If ymeas is used as the feedback signal to Ktheta (or in the commented out line of Kstate), then Ktheta*P(s)/s won't show up in the numerator of H(s), but the Kvel*P(s) term will still dominate."

Hmmm, I'll have to dig into this more. Since I see FB actively changing its signal to mirror and cancel FF, wouldn't that continue regardless of Gff (F)? That still feels like LQI trying to maintain its own disturbance-rejection profile -- but that theory might be wrong.

[Paul]

Yes, if the state variable theta is defined as the angle of the platform relative to one reference, like the local level, and the feedback to the integrator is the angle of the platform relative to the target (which is the ouput disturbance) that should be regulated to zero, then the Simulink model is properly capturing this situation as far as I can tell.

As for that last paragraph, I'm not sure what you mean by "mirror and cancel." In either case, the input to the

Repeating the code

load FF_check
plant2.C = [plant2.C;eye(2)];
plant2.OutputName = {'y' 'x1' 'x2'};
gff2.InputName = 'd';
gff2.OutputName = 'uffin';
sensor = sumblk('ymeas = y + d');
Ki     = ss(K_lqi2(3),'InputName','ei','OutputName','ui');
Kstate = ss(K_lqi2(1:2),'InputName',{'x1' 'x2'},'Outputname','ustate');
Kfb    = tf(realp('Kfb',1),'InputName','ufbin','OutputName','ufb');
Kff    = tf(realp('Kff',1),'InputName','uffin','OutputName','uff');
sys = connect(...
    sumblk('e = r - ymeas'), ...
    tf(1,[1 0],'InputName','e','OutputName','ei'), ...
    Ki, ...
    Kstate, ...
    plant2, ...
    sensor, ...
    gff2, ...
    sumblk('ufbin = -ui - ustate'), ...
    sumblk('u     = uff + ufb'), ...
    Kfb, ...
    Kff, ...
    {'r' 'd'}, ...
    {'y' 'ymeas' 'u' 'uff' 'ufb'});

Beause the output ymeas is essentially the same for Kff on or off (as justified above), we should expect the input to the plant to be essentially the same as well.

Plot the u responses for both cases and we see that the input to the plant is essentially the same for both cases, except for the very initial transient due to the direct feedthrough of gff2.

[uon,ton]   = step(sampleBlock(sys({'u' 'uff' 'ufb'},'d'),'Kfb',1,'Kff',1));
[uoff,toff] = step(sampleBlock(sys({'u' 'uff' 'ufb'},'d'),'Kfb',1,'Kff',0),ton);
figure
plot(ton,uon(:,1),toff,uoff(:,1)),grid
xlim([0 .5e-4])

Now if the input to the plant is (essentially) the same for both cases, we'd definitely expect ufb to be different for both cases. In one case we have u = ufb and in the other we have u = ufb + uff.

As far as what gff2 should be, consider

P      = zpk(plant2('x2','u'));
Kvel   = K_lqi2(2);
Ktheta = K_lqi2(1);
s      = zpk('s');
% H = (F*P/s + 1 + Kvel*P + Ktheta*P/s)/(1 + Ki*P/s^2 + Ktheta*P/s + Kvel*P);

To make H(s) = 0, we'd have mathematically:

F = -s/P*(1 + Kvel*P + Ktheta*P/s);

Or more cleanly

F = -minreal(s/P) - s*Kvel - Ktheta

F =
 
  From input "x2" to output:
  -7.0362e-06 (s+4.155e06) (s+1.713e06)
 
Continuous-time zero/pole/gain model.

So just make gff2 the low frequency gain to kill the effect of the disturbance on ymeas (at least mathematically, the plant input might be very unrealistic).

gff2.Num = freqresp(F,0);
gff2.Den = 1;

sys = connect(...
    sumblk('e = r - ymeas'), ...
    tf(1,[1 0],'InputName','e','OutputName','ei'), ...
    Ki, ...
    Kstate, ...
    plant2, ...
    sensor, ...
    gff2, ...
    sumblk('ufbin = -ui - ustate'), ...
    sumblk('u     = uff + ufb'), ...
    Kfb, ...
    Kff, ...
    {'r' 'd'}, ...
    {'y' 'ymeas' 'u' 'uff' 'ufb'});
figure
step( ...
    sampleBlock(sys({'ymeas'},'d'),'Kfb',1,'Kff',1), ...
    sampleBlock(sys({'ymeas'},'d'),'Kfb',1,'Kff',0));
legend('Kff = 1','Kff = 0')

[John]

Thanks Paul.

1) "As for that last paragraph, I'm not sure what you mean by "mirror and cancel." In either case, the input to the "

(curious how that sentence was going to end)

Here's what i meant. For the plot of u_FB from above, it mirrors u_FF to cancel out FF's contribution. Your u commands didnt show this overshoot command in u_FB, and only had the initial negative-going peak. I'm not sure why they're different, since we're implementing the same thing...

2) "It seems to me that the strategy for selecting F(s) was to ensure that F(s)*P(s)/s + 1 is approximately zero at low frequency. However, that only takes care of two terms in the numerator of H(s)."

If the realizable plant pseudoinverse is used to define F (F ~= -1 / P), doesn't that mean F*P = ~1 for all frequencies, not just near DC?

Note: I could have removed the integrator following F, and then also the derivative following the disturbace, but added them there to be consistent with the physical measurement. But they cancel out, so it seems it'd be F(s)*P(s) = 1 for all frequencies. But yes, it should settle to -1 ultimately (canceling out d), which is perhaps all that you were saying by this.

3) "However, the numerator of H(s) is dominated by Kvel*P(s) and Ktheta*P(s)/s"

"The closed-loop transfer function from d to ymeas is then:

H = (F*P/s + 1 + Kvel*P + Ktheta*P/s)/(1 + Ki*P/s^2 + Ktheta*P/s + Kvel*P);"

To understand more, how did you find this H(s)? I had something else:

syms P F Kw Kth Ki s;
% uff is the command, passing through integrator: ff = F/s
% Integrator added if using s*dist to get dist angvel
syms uff fbGainPath;
uff = F/s;
innerfbGainPath = Kw*P*s + Kth*P;
innerFwdPath = P*(1+uff);


% SP --> y
% Inner loop: fwd path is plant and FF input, and fb path is gains*states
Hin = innerFwdPath / (1 + innerFwdPath*innerfbGainPath);
% Outer loop: fwd path is integrator path, and inner loop. Fb path is 1.
Hout = Ki/s * Hin / (1 + Ki/s * Hin);
simplify(Hout)

ans = 

% d --> y
% Inner loop: fwd path is 1; fb path is gains*states, plant, and FF input
innerFwdPath = 1;
innerfbGainPath = (Kw*P*s + Kth*P) * P*(1+uff);
Hin = innerFwdPath / (1 + innerFwdPath*innerfbGainPath);
% Outer loop: fwd path is 1; fb path is integrator path + and inner loop.
outerFwdPath = 1;
outerFbPath = Ki/s * Hin;
Hout = outerFwdPath / (1 + outerFwdPath*outerFbPath);
simplify(Hout)

ans = 

I'll assume your TF is correct: I do see that that the TF is dominated by gains*states, but in my SL model i also see FB actively canceling out a large FF command contribution, seemingly in order to maintain its designed rejection trajectory.

If H(s) was dominated by those terms, then why would FB change its command (when you run my model) to cancel out FF? FB would -- i think -- simply not need to change u_FB, because u_FF contributes so little compared with dominant terms. Instead, FB actively cancel FF's relatively large contribution, at least in my SL model.

Is this not because LQI has a designed rejection trajectory, and u_FF looks like a disturbance?

4) "Plot the u responses for both cases and we see that the input to the plant is essentially the same for both cases, except for the very initial transient due to the direct feedthrough of gff2."

Odd...that's not what the simulink model's output shows... However, they should be the same given how you closed the loop using connect().

What do you see for the command signals when you run SL model?

There have been times when i've posted here in some confusion, and it turned out the that ML/SL had a bug that caused some issues -- so double checking what you see.

5) "So, in the current form of the system (d only affecting feedback to integrator), it seems like F(s) should be selected such that it is realizable and satisfies

F(s)*P(s)/s + 1 + Kvel*P(s) + Ktheta*P(s)/s = 0 at low frequency."

"So just make gff2 the low frequency gain to kill the effect of the disturbance on ymeas"

How come this isn't F(s)*P(s)/s + 1 + Kvel*P(s)*s + Ktheta*P(s) = 0?

In any case: I implemented F = -minreal(s/P) - s*Kvel - Kthetain SL to check (since intuitively i don't follow why the additional terms are needed), but wasn't seeing the same results as your step response plot.

The SL plot below shows FB+FF, and I still see FB canceling out FF's command; the response is a bit less damped than FB-only. I added a pair of fast poles to make it implementable:

-7.0362e-06 * (s+4.155e06) * (s+1.713e06) * 5000^2/(s+5000)^2

If i slow down the added poles (to say 1000 rad/s), the response is exactly the same as FB-only.

[Paul]

1) Sorry about the orphaned sentence on the "mirror and cancel" thing. I was headed down the path to show that the plant input, u, is essentially the same with Kff = 1 and Kff = 0. I guess I got ahead of myself to show the result and then didn't come back to type the text.

It does appear that we have different results. I'm not running your Simulink model. I hope my Matlab code is consistent with what you've modeled in Simulink. I do see that it looks like your model is exciting the system such that "y after disturbance" has a step jump to a value of 100 very shortly after t = 0, without any change in uff or ufb that time, and then something else happens at t = 0.0006 to drive "y after disturbance back to zero." Have you tried running your model as I've done here, i.e., starting with all states at zero and then apply a unit step input disturbance?

2) But as I showed, making F(s) = -1/P(s) isn't the right approach to mitigated the effect of d on y, at least based on my understanding of the system.

3) ufb will do whatever it needs in accordance with the closed-loop dynamics. With gff2 as defined, the ufb with Kff = 1, let's call it ufbwithff, will be whatever it nees to be such that ufbwithff + uff = ufbwithoutff (roughly).

4) I'm not running your SL model. If it will help, I might be able to develop a SL model of what I think the system is and compare it to the ML code

5) The only thing to implement is to make gff2 a constant value

[John]

1) "and then something else happens at t = 0.0006 to drive "y after disturbance back to zero."

That's a feedback sensor delay in the SL model: it's in the feedback path. Also, you can see a delay in creating the FF command: that's the ff sensor delay. The delays can be removed without significant change in behavior.

"Have you tried running your model as I've done here, i.e., starting with all states at zero and then apply a unit step input disturbance?"

Yes -- actually, I thought that's what i've been running. I have changed it to step 1 instead of step 100, but that only scales everything by 100 :) The SL plots i've been showing have been all with steps of 100, starting from 0-state.

2) Also, check out the new plot i added to (5) in my prev post, which shows the FB command again canceling the FF command. The plot looks the same despite a constant gff, ie FB cancels FF.

3) "The only thing to implement is to make gff2 a constant value"

How come? It looked like you had F = -7.0362e-06 (s+4.155e06) (s+1.713e06)

But perhaps I misunderstood what you were suggesting.

"With gff2 as defined, the ufb with Kff = 1, let's call it ufbwithff, will be whatever it nees to be such that ufbwithff + uff = ufbwithoutff (roughly)."

Agreed, and that seems to be the issue...in general, inverse plant dynamics is the right thing for FF (as evidenced by the nice response with FF-only).

Even if gff is just a constant value, the FB command cancels the FF command as before; see below. It doesn't seem to matter what gff is, since it seems LQI tries to cancel it out to maintain its dynamics.

4) " But as I showed, making F(s) = -1/P(s) isn't the right approach to mitigated the effect of d on y, at least based on my understanding of the system."

"I hope my Matlab code is consistent with what you've modeled in Simulink."

I'm not sure that we have the same system: i looked at the loop TF and seem to get something different (sorry -- I added that part later in a post edit of my previous post).

Also, my u_FF and u_FB don't seem to match what you have, even when i remove any Group Delays.

"I'm not running your SL model. If it will help, I might be able to develop a SL model of what I think the system is and compare it to the ML code"

I'd be curious what you see if you run the model i attached, to see if the response is consistent with what you expect given your ML code. (Comment-through the two Group Delays to match your code). The model is fairly basic and captures what we discussed, so I'm not sure if it's worth building up another one...

[Paul]

1) why not get rid of all of the delays so we're comparing the same models.

3) With my definitions (P(s) is the transfer fucntion from u to ang_vel as I stated previously), we'd ideally want the feedforward compenation, F(s), to be such that

(F*P/s + 1 + Kvel*P + Ktheta*P/s) = 0

which is exactly the F that I solved for:

F =

From input "x2" to output:

-7.0362e-06 (s+4.155e06) (s+1.713e06)

I approximated this as F(s) = -7.0362e-6 * 4.155e6 * 1.713e6 (to full precision in each term) by evaluating freqresp(F,0)

I derived the transfer function H(s) from this diagram, which is the system as I understand it (w/o delays)

There are two forward paths from d to ym.

The forward path gain of the first is: F1(s) = gff2(s)*P(s)/s

The forward path gain of the second is: F2(s) = 1*(1 + Kvel*P(s) + Ktheta*P(s)/s)

Note that for F2 we have to account for the fact that the path from d to ym does not break the two inner loops.

The characteristic equation is just from the three loops: D(s) = 1 + Ki*P(s)/s^2 + Ktheta*P(s)/s + Kq*P(s)

[John]

I see, maybe that's the difference: P is a command-to-position plant (not velocity). Imagine that P is a spring, so a command input causes a certain change (but not indefite motion). The difference is that gff acting on P would immediately change y, vs needing to integrate P to get y.

That might also explain the H(s) differences...

[Paul]

Here's the plant that you provided

load FF_check
plant2

plant2 =
 
  A = 
                theta      angVel
   theta            0           1
   angVel  -1.421e+05        -240
 
  B = 
                   u
   theta           0
   angVel  1.421e+05
 
  C = 
           theta  angVel
   theta       1       0
 
  D = 
          u
   theta  0
 
Continuous-time state-space model.

As we can see, D = 0 so there can't ever be an immediate change in theta in response to a change in u, as can be seen from the step response

figure
step(plant2)

Here is gff2

gff2

gff2 =
 
  -3.198e09 s^2 - 7.675e11 s - 4.545e14
  -------------------------------------
  1.421e05 s^2 + 1.607e10 s + 4.545e14
 
Continuous-time transfer function.

It is a proper transfer function, so placing it in series with the plant also can't cause an immediate change in theta in response to a step input into gff2

figure
step(plant2*gff2)

In my analysis and in the block diagram I showed, P(s) is NOT tf(plant2). Rather, it is the transfer function from U(s) to AngVel(s), as shown in the diagram, stated previously, and computed as such in this comment. thetadot = angVel, as shown in plant2 and represented in my analysis and block diagram.

In your analysis, it looks like P(s) is intended to be the transfer function from U(s) to Theta(s). But that's ok. The results should be the same as long as the analysis is done properly. Feel free to provide a complete block diagram of the system, even handwritten, if you still think that my representation (based on the portion of the Simulink model shown in an upstream comment) is incorrect.

Here's the block diagram using what I think is your definition of P(s), which I'll call P2(s), and is the transfer functiom from U(s) to Theta(s), i.e., P2(s) = tf(plant2)

Based on this diagram, the transfer function H2(s) = Ymeas(s)/D(s) is:

H2(s) = ( F(s)*P2(s) + 1*(1 + Kvel*P2(s)*s + Ktheta*P2(s) ) / (1 + Ki*P2(s)/s + Kvel*P2(s)*s + Ktheta*P2(s)

which is the same as what I showed above

H = (F*P/s + 1 + Kvel*P + Ktheta*P/s)/(1 + Ki*P/s^2 + Ktheta*P/s + Kvel*P);

i.e., H(s) == H2(s) because P2(s) = Theta(s)/U(s) = AngVel(s)/U(s)/s = P(s)/s.

[John]

1) I have more to ask later :), but was hoping to clarify Matlab nomenclature first, since that may change my questions...

"Define the controller as u = kfb*(r - y) - B(2)*wmeas, i.e., augment the feedback control with the negative of the disturbance input to the plant."

I see that you're using LQR by way of example (since this post is about LQI), but I thought Matlab's LQR uses the inner loop here where states are in a direct feedback (ignoring integrator part of course) and without reference error, or here showing an LQR loop. To understand more, why are you using K*(r-y), instead of K*y? Or does matlab have different LQR configurations that are tracking vs regulating?

2) "In my analysis and in the block diagram I showed, P(s) is NOT tf(plant2). Rather, it is the transfer function from U(s) to AngVel(s), as shown in the diagram, stated previously, and computed as such in this comment. thetadot = angVel, as shown in plant2 and represented in my analysis and block diagram.

In your analysis, it looks like P(s) is intended to be the transfer function from U(s) to Theta(s)."

"B =

u

theta 0

angVel 1.421e+05"

Since the output Y is theta as given by C in Matlab's formulation, how come this isn't considered u --> theta in Matlab's formulation?

My block diagram used the TF as u-->theta (Y(s) / U(s) = pos / cmd) in the .slx file i attached, so the states output x used by LQI in my .slx are x1=theta and x2=angVel.

I'm trying to understand why the B above would presume u --> angvel TF (as you noted), vs what i thought would cause u-->theta: the input stimulating velocity but leading to a position change (eg a spring that settles).

[Paul]

As for (1), my first comment in this thread wasn't intended to have anything to do with LQR or LQI. All I was trying to show was that a disturbance at the input to the plant could be exactly cancelled in the control input. It wasn't until later in the thread that it became clear to me that the disturbance is added at the output of the plant.

As for (2), I've read your comment a few time and I'm still not sure what you're asking.

Throughout, I've been using the variable plant2 that you uploaded

load FF_check
plant2

plant2 =
 
  A = 
                theta      angVel
   theta            0           1
   angVel  -1.421e+05        -240
 
  B = 
                   u
   theta           0
   angVel  1.421e+05
 
  C = 
           theta  angVel
   theta       1       0
 
  D = 
          u
   theta  0
 
Continuous-time state-space model.

As you noted, B = [0 ; 1.421e5]. Is that not what you expected it to be?

The state variables are theta and angVel, and those are the states used for the state feedback portion of the control.

For purposes of deriving the closed loop transfer function from d to ymeas ( = d + theta), the transfer function form of plant2 is

P2 = tf(plant2)

P2 =
 
  From input "u" to output "theta":
         1.421e05
  ----------------------
  s^2 + 240 s + 1.421e05
 
Continuous-time transfer function.

and it is this P2(s) that would be used in the block diagram of this comment to derive the closed-loop transfer function H(s) = Ymeas(s)/D(s). Note that in that block diagram the feedback variables to the state feedback gains are theta and s*theta = angVel, as they should be.

In the block diagram in this comment I used a plant model P(s) that is the transfer function from u to angVel and then used an integrator to get from angVel to theta. That transfer function would be

P = tf(ss(plant2.A,plant2.B,[0 1],0))  % A and B don't change!

P =
 
        1.421e05 s
  ----------------------
  s^2 + 240 s + 1.421e05
 
Continuous-time transfer function.

which, of course, satisfies P(s)/s = P2(s). In this block diagram, the feedback variables to the state feedback gains are still theta and angVel, as they should be.

The two block diagrams are equivalent and yield the exact same H(s).

[John]

Thanks Paul. My question was confusing because i grammatically misunderstood something you wrote :) Your latest explanation said the same thing a different way, so it clarified and makes sense now.

1) Does s = zpk('s') that's in your code, differ in usage/capabilities from s = tf('s') that i typically use?

2) "Feel free to provide a complete block diagram of the system, even handwritten..."

Your updated diagram matches the block diagram in my attached .slx file :)

Mine has switches, and uses states x1 and x2 to get vel as opposed to your P2*s; but functionally they match. I also removed the delays from y and x to match yours:

3) What is this nomenclature doing, that you had above?

P = plant2('x2','u')

I'm not familiar with it, and don't know what to search for. When I tried, there was an error (Error using (): The name "x2" does not match any channel or group.)

I repeated with

P = plant2('theta','u');

which also gave the error, though I see 'theta" as one of the SS name params.

3) "Based on this diagram, the transfer function H2(s) = Ymeas(s)/D(s) is: <...> which is the same as what I showed above"

Just to lean more: did you find H(s) using matlab and connect(), or another method?

The script i posted in this comment I thought captured compute F2(s): first, the inner loops' TF including Ktheta*s and Kvel, then the outer loop including Ki/s.

But, your F2(s) path, which includes the outer-most loop with Ki, doesn't seem to include Ki/s, so I'm not following F2(s) = 1*(1 + Kvel*P2(s)*s + Ktheta*P2(s)

How come F2(s) doesn't include -Ki/s * P2(s)?

4) Assuming your TF is correct (i'm assuming mine isn't), I tried your proposals for H(s) --> 0:

• gff2=const, ie F evaluated at 0 ("The only thing to implement is to make gff2 a constant value"...ie -5.008e7), and
• gff2=function of Ktheta and Kvel: "F = -7.0362e-06 (s+4.155e06) (s+1.713e06)", but with an implementable version with 2x fast poles: *(1e8^2 / (s+1e8)^2), with ~same result as the constant

The result is as you said: FB+FF now has a faster output as FB-only.

But there are unexpected outcomes:

(a) A FF-only of either method above is ~300x slower than FF-only of -1/P (my original method).

(b) With FB+FF, FB cmd still tries to mostly cancel the u_FF, ie at steady-state, FB is ~ -5.008e7, and FF is ~ +5.008e7.

In an ideal, properly designed FB+FF system, FF does the bulk of the work, and FF would just be the residual error trim. Ie, FB+FF is faster than FF-only (which does occur in this case), but also FF-only works well by itself without the FB trim.

Instead, in this case FF does not work well without the feedback trim.

Also, with any delays in the system, the stability of this system seems heavily impacted since FF needs to be more precisely computed given knowledge of FB. This is also different than what i've typically encountered: FF is designed independently from FF, then added to FB, with typically excellent FB+FF results.

(c) At steady-state, FB+FF = 203, which seems concerning stability-wise with practical implementation: precision needs to be 1/1000 % (200 out of 5e7).

You also mentioned practical considerations: "So just make gff2 the low frequency gain to kill the effect of the disturbance on ymeas (at least mathematically, the plant input might be very unrealistic)."

Overall, since typically FF is designed independently from FB and they are usually just tied together without general loss of stability, I'm still unclear why LQR/LQI are different (than say PID). It still seems like there's something about the LQI loop structure (full-state FB, vs PID that's all error-based FB) or paramaterization that makes this less straightforward.

[Paul]

@John, your last two posts were posted as new answers. I hope you don't mind me moving them to comments in this single answer thread.

1) zpk('s') just uses the zpk form instead of the tf form. See zpk for more info if you're not familiar. AFAIK, usage is the same as tf. Generally speaking, ss form is preferred for computations, though zpk and ss are sometimes better for visual display. I think also that zpk is preferred to tf for computations if having to choose between one or the other. Using the Right Model Representation

2) So I guess you're results should match mine as far as 'y after dist' in response to 'dist_pos', which is what I called the ymeas in response to d.

3) We can use output names and input names to select specific parts of MIMO models. In the code in the comment where I did that I had added another output to the plant and gave that output an outputname of 'x2'. For example, suppose we have a ss system with two input and two outputs

sys = zpk(rss(1,2,2));sys.InputName = {'u1' 'u2'}; sys.OutputName = {'y1' 'y2'};

If we want the system from input1 to output2 then we can index numerically

sys(2,1)

ans =
 
  From input "u1" to output "y2":
  -0.3682
 
Static gain.

Or

sys('y2','u1')

ans =
 
  From input "u1" to output "y2":
  -0.3682
 
Static gain.

Use cell arrays for more than one input and/or output. Here, we can reverse the inputs and outputs

sys({'y2' 'y1'},{'u2' 'u1'})

ans =
 
  From input "u2" to output...
        0.033242
   y2:  ---------
        (s+3.291)
 
   y1:  0
 
  From input "u1" to output...
   y2:  -0.3682
 
   y1:  0.42888
 
Continuous-time zero/pole/gain model.

Can just as easily use numeric indexing, but the char indexing is more clear IMO.

I'm SHOCKED that char indexing into MIMO models is not documented here Select Input/Output Pairs in MIMO Models, but it is used here Robust Control of Active Suspension.

3) I found H(s) by hand based on the block diagram.

The block diagram has three loops all touching each other. There are two direct paths from d to ymeas. Recall from Mason's Rule that the forward path gain is the direct path gain multiplied by the characteristic equation after removing all loops touched by the direct path. Path1 touches all three loops, so the forward path gain is just the direct path gain F(s)*P2(s). Path2 has a direct path gain of 1, but the direct path doesn't touch the two inner loops, hence the path gain for Path2 is 1*(1 + Kvel*s*P2(s) + Ktheta*P2(s)).

4) "The result is as you said: FB+FF now has a faster output as FB-only. " I hope you're seeing in your Simulink model that FB+FF is much, much faster as I showed above.

4a). In the absence of feedback, the transfer function from d to ymeas is simply:

Ymeas(s)/D(s) = 1 + Gff2(s)*P(s)

where P(s) = Y(s)/U(s).

If you make Gff2(s) = -1/P(s)*SuperFastPoles(s), that will be a much faster response then 1 + K*P(s), where K is constant gain (-5.008e7).

4b). See comments on 4c for response to the first part.

As to the stability of the system, I would imagine there would be sensitivity to small delays in a system with such high loop bandwidth (which is not impacted by the feedforward control), but I'm not sure what you mean by stability in the sense of "FF needs to be more precisely computed given knowledge of FB."

4c) My practical considerations are more focused on the transient input into the plant. I suspect that the feedforward term is going to cause an extremely large magnitude input with a very high rate and the physical sytem might not be able to respond to such an input in the way the mathematical model predicts.

load FF_check
% Use the feedforward compensatin derived above
gff2.Num = -7.0362e-06 * (4.155e06) * (1.713e06);
gff2.Den = 1;
plant2.C = [plant2.C;eye(2)];
plant2.OutputName = {'y' 'x1' 'x2'};
gff2.InputName = 'd';
gff2.OutputName = 'uffin';
sensor = sumblk('ymeas = y + d');
Ki     = ss(K_lqi2(3),'InputName','ei','OutputName','ui');
Kstate = ss(K_lqi2(1:2),'InputName',{'x1' 'x2'},'Outputname','ustate');
Kfb    = tf(realp('Kfb',1),'InputName','ufbin','OutputName','ufb');
Kff    = tf(realp('Kff',1),'InputName','uffin','OutputName','uff');
sys1 = connect(...
    sumblk('e = r - ymeas'), ...
    tf(1,[1 0],'InputName','e','OutputName','ei'), ...
    Ki, ...
    Kstate, ...
    plant2, ...
    sensor, ...
    gff2, ...
    sumblk('ufbin = -ui - ustate'), ...
    sumblk('u     = uff + ufb'), ...
    Kfb, ...
    Kff, ...
    {'r' 'd'}, ...
    {'y' 'ymeas' 'u' 'uff' 'ufb'});

The feedforward term significantly speeds up the response of the output

% step response from disturbance to ymeas
figure
step(sampleBlock(sys1('ymeas','d'),'Kfb',1,'Kff',1), ...
     sampleBlock(sys1('ymeas','d'),'Kfb',1,'Kff',0))

at the expense of driving the system with a very large input. Can the physical system tolerate that input?

% step response from disturbance to u
figure
step(sampleBlock(sys1('u','d'),'Kfb',1,'Kff',1), ...
     sampleBlock(sys1('u','d'),'Kfb',1,'Kff',0))

How did you get that the steady state value of the input to the plant is uss = 203? Based on the definiton of plant2, the feedback control forces the steady state input to the plant, uss, in response to a unit step input disturbance, d=1, to be uss = -1. This result follows from the plant2 model equation for the angular acceleration

thetadotdot = A21*theta + A22*thetadot + B21*u

At steady state, in response to a step disturbance, we have theta = -d, thetadot = 0, and thetadotdot = 0, which means

0 = -A21*d + B21*uss -> uss = A21/B21*d = -1

for d = 1 and the given values of A21 and B21

plant2.A(2,1)/plant2.B(2,1)

ans = -1

We can verify this result by computing the final value of the plant input from a unit step disturbance with FF on and off

[freqresp(sampleBlock(sys1('u','d'),'Kfb',1,'Kff',1),0);
 freqresp(sampleBlock(sys1('u','d'),'Kfb',1,'Kff',0),0)]

ans = 2×1

-1.0000 -1.0000

The equation for the input to the plant is

u = GFF2*d - Kvel*thetadot - Ktheta*theta + Ki*ei

The steady state input for step disturbance with integral control must satisfy

uss = GFF2*d - Kvel*0 + Ktheta*d + Ki*ei -> uss = (GFF2 + Ktheta)*d + Ki*ei

But, we know that uss = -1, and if d is unit step

-1 = GFF2 + Ktheta + Ki*ei

where the first term on the RHS if what you're calling FF and the sum of the second two terms is what you're calling FB. Therefore, in steady state

-1 = FF + FB

In the absence of feedforward control, the steady state feedback would be

FB = -1

With the feedforward control, we'd have

FB = -1 - FF

which is essentially FB = -FF as you've observed because FF is enormous compared to unity for the choice of Gff2.

All of the results make sense (at least to me) in the context of the control system as implemented.

For this problem, for one thing, the state feedback configuration is different than PID insofar as with PID there's only one loop. So, with PID implemented in K(s), the closed loop transfer function from D(s) to Ymeas(s) would be

H(s) = (1 + Gff2(s)*P(s))/(1 + K(s)*P(s))

Now, setting Gff2(s) = -1/P(s)*SuperFastPoles(s) would do what I think you were expecting at the outset.

[John]

"your last two posts were posted as new answers. I hope you don't mind me moving them to comments in this single answer thread."

Apologies -- and thanks for fixing :)

Thanks to your explanations this is all making much more sense:

"I think also that zpk is preferred to tf for computations if having to choose between one or the other. Using the Right Model Representation"

Thanks for the explanation -- and that's a great link, that i'm surprised i haven't run across at some point. I'll convert my libraries to zpk().

"So I guess you're results should match mine as far as 'y after dist' in response to 'dist_pos', which is what I called the ymeas in response to d."

"The result is as you said: FB+FF now has a faster output as FB-only. " I hope you're seeing in your Simulink model that FB+FF is much, much faster as I showed above."

Yes, they seem to match exactly :) [when F(s) input effort is unconstrained]

"We can use output names and input names to select specific parts of MIMO models. <...> Can just as easily use numeric indexing, but the char indexing is more clear IMO. <...> I'm SHOCKED that char indexing into MIMO models is not documented here Select Input/Output Pairs in MIMO Models"

That's very useful. I'm familiar with numeric indexing but I'm not sure why the label indexing isn't in the basic state space documention.

"In the absence of feedback <...>, Ymeas(s)/D(s) = 1 + Gff2(s)*P(s)

If you make Gff2(s) = -1/P(s)*SuperFastPoles(s), that will be a much faster response then 1 + K*P(s), where K is constant gain (-5.008e7)."

"For this problem, for one thing, the state feedback configuration is different than PID insofar as with PID there's only one loop. So, with PID implemented in K(s), the closed loop transfer function from D(s) to Ymeas(s) would be <...>

Now, setting Gff2(s) = -1/P(s)*SuperFastPoles(s) would do what I think you were expecting at the outset."

Ah; I'm following now and this all makes sense -- thanks. To restate in my own words:

For full-state FB, while keeping FF On at all times: if turning FB On or Off for any reason, F(s) needs to dynamically change to match the FB loop characteristics such that Ymeas(s)/D(s) = 0 (ie smoothly switching between -1/P and -1/P +...). Adding FF in PID is straightforward because of where the FF cmd is injected relative to the error-based feedback (in most constructions has no inner loops).

So, the general case is: FF needs to take loop structure into account, and people are simplifying (and assuming error-based feedback) when saying that FB and FF can be designed independently with the expectation of performance increase after connecting both. This isn't (in general) true, as seen in this post.

"I found H(s) by hand based on the block diagram. <...> There are two direct paths from d to ymeas. Recall from Mason's Rule that the forward path gain is the direct path gain multiplied by the characteristic equation after removing all loops touched by the direct path. <...>"

Ah, okay. Block reduction (which was the script I posted) should reach the same result, but i must have made an error. So I just went through the full example using Mason's (not that i didn't trust your result...for my own understanding) and also reached the same result of F(s) = -1/P + Ktheta - Kvel*s, to make H(s) -->0. One question though:

1) Why is it sufficient to say H(s) = 0, when there's no performance criteria attached? It seems this means converge to 0 for all frequencies as t --> inf, but it could be quick or slow evolution.

2) "The steady state input <...> must satisfy

uss = GFF2*d - Kvel*0 + Ktheta*d + Ki*ei -> uss = (GFF2 + Ktheta)*d + Ki*ei

In the absence of feedforward control, the steady state feedback would be FB = -1.

With the feedforward control, we'd have FB = -1 - FF

which is essentially FB = -FF as you've observed because FF is enormous compared to unity for the choice of Gff2."

This makes sense...but this is all in DC, no? What i see is that FB ~= -FF even for high initial impulses (per my images in previous posts), ie it seems this holds for all transient time, ie thetadot != 0.

3) "We can verify this result by computing the final value of the plant input from a unit step disturbance with FF on and off

[freqresp(sampleBlock(sys1('u','d'),'Kfb',1,'Kff',1),0);

freqresp(sampleBlock(sys1('u','d'),'Kfb',1,'Kff',0),0)]"

I follow your reasoning, but to check, why use freqresp()? Would dcgain(sampleBlock(sys1(...))) not accomplish the same thing?

[Paul]

1) From the outset, I was under the impression based on your definition of Gff2(s) was that you were trying to define Gff2(s) so that, ideally, the numerator of H(s) = Ymeas(s)/D(s) = 0, but for practical reasons you included two very fast poles in Gff2(s). As we've discovered, the original definition of Gff2(s) did not take into account all of the paths from the disturbance, d, to ymeas. I just followed what I thought the design goal is, which is to drive the ymeas to zero as quickly as possible. Now that you have the correct expression for H(s), you can modify Gff2(s) to change the response to try to meet other performance criteria.

2) I don't observe the same behavior.

load FF_check
% Use the feedforward compensation derived above
gff2.Num = -7.0362e-06 * (4.155e06) * (1.713e06);
gff2.Den = 1;
plant2.C = [plant2.C;eye(2)];
plant2.OutputName = {'y' 'x1' 'x2'};
gff2.InputName = 'd';
gff2.OutputName = 'uffin';
sensor = sumblk('ymeas = y + d');
Ki     = ss(K_lqi2(3),'InputName','ei','OutputName','ui');
Kstate = ss(K_lqi2(1:2),'InputName',{'x1' 'x2'},'Outputname','ustate');
Kfb    = tf(realp('Kfb',1),'InputName','ufbin','OutputName','ufb');
Kff    = tf(realp('Kff',1),'InputName','uffin','OutputName','uff');
sys1 = connect(...
    sumblk('e = r - ymeas'), ...
    tf(1,[1 0],'InputName','e','OutputName','ei'), ...
    Ki, ...
    Kstate, ...
    plant2, ...
    sensor, ...
    gff2, ...
    sumblk('ufbin = -ui - ustate'), ...
    sumblk('u     = uff + ufb'), ...
    Kfb, ...
    Kff, ...
    {'r' 'd'}, ...
    {'y' 'ymeas' 'u' 'uff' 'ufb'});

System from d to u, uff, and ufb.

sys = sampleBlock(sys1({'u','uff','ufb'},'d'),'Kfb',1,'Kff',1);

Verify DC gain

format long
dc = dcgain(sys);

As expected from above, the steady state value of plant input is -1.

dc(1)

ans = 
    -1

As expected, at steady state the sum of uff and ufb is -1, equal to u.

dc(2)+dc(3)

ans = 
    -1

Step response from d to u, uff, and ufb.

figure
step(sys)

The middle plot shows that uff is a constant as should be the case for a step input disturbance and Gff2 being a constant. The bottom plot shows that ufb has a dynamic transient response. I don't see how it would be the case that ufb(t) is approximately -uff(t), except as ufb approaches steady state.

3) Yes, dcgain would accomplish the same thing. I just forgot about that function.

[John]

Thanks Paul.

(Post edited as i realized I was wrong on (1) )

1) In the numerator, i had (F*P + <...>) instead of F*P/s. <-- Edit: you also had this. NVM.

On a related note, it's interesting how the full F(s) doesn't reach the performance possible with just the constant F(0) = -5e7, since the fast poles seem to slow it down. Or if the poles are made faster to reduce their impact, there is significant overshoot (eg above 3e6 for the fast pole locations). I'm not sure why the full-order F model is worse than the const version -- at least at the limit.

When not at the limit, the extra fast poles of F can be tuned to adjust the sytem settling performance. I'm not sure what other method is available, though, beyond that.

2) "I don't see how it would be the case that ufb(t) is approximately -uff(t), except as ufb approaches steady state."

Agreed for F(0) = constant, but check out the non-const version of Gff2 * fastPoles. Eg:

-7.0362e-06 * (s+4.155e06) * (s+1.713e06) * (1e6^2/ (s+1e6)^2)

This will show a symmetric response when FF and FB are active:

And if FF is Off over the same time interval, ie FB only, shows what FB normally would be without FF. So FB does seem to significantly change shape based on FF, which you did show in your last post when you showed FB = -1 - FF, which these plots support as t-->inf. So, only posting this item (2), to confirm that i also see the limit u = -1 (not 200 as i originally said -- I just didn't look far enough into the future)

[Paul]

I'm glad we seem to have it all straightened out.

[John]

Thanks again for your help -- and (as always) the additional Matlab methods you show to analyze such problems, are quite helpful.

[John]

Hi @Paul, just going back to try and understand something about the step response of H(s) = D(s) / Ym(s).

In your post here (if link doesn't work, it's 27 Dec 2023 at 22:37), you plotted the step response d --> y_measured, after using connect() to create the FB+FF system.

In another post here (28 Dec 2023 at 13:47), you showed H(s) = Ym(s) / D(s), ie the same thing as the first post i linked:

"Based on this diagram, the transfer function H2(s) = Ymeas(s)/D(s) is:

H(s) = ( F(s)*P(s) + 1*(1 + Kvel*P(s)*s + Ktheta*P(s) ) / (1 + Ki*P(s)/s + Kvel*P(s)*s + Ktheta*P(s)"

(fyi i changed eg P2,H2 to P,H, just for clarity in this post)

Later i'd also verified H(s) = Ym(s)/D(s), as well as F(s) = -1/P - Ktheta - Kvel*s, so i'll assume they're correct.

It seems your connect() method's step response should yield the same result as step(H), -- is that right?

I ask because step(H) seems to not converge, whereas your connect()'s d-->ym plot clearly does converge to -1, so I'm trying to understand if the connect() block-diagram method should yield the same result as the freq-domain H(s), which i would have thought they should.

[Paul]

The block diagram in the second linked comment form the input to the plant as:

u = ui - ufb

whereas the block diagram represented by connect in the first linked comment forms the input to the plant as

u = -ui - ufb

So the algebraic approach to forming H has to use the negative of Ki and compared to connect.

With connect (first linked comment) I was trying to be faithful to the block diagram in your Simulink digram, which is exactly how the system is shown on lqi.

With algebraic approach, I was thinking of the problem from a typical control system how integral control on an error signal.

Accounting for that difference, the two approaches match.

I should have been more clear about that. If I get a chance I'll go back and edit the relevant comments and add that clarification.

If that doesn't solve the problem, I guess you'll have to post your code to see what's happening.

load FF_check
plant2.C = [plant2.C;eye(2)];
plant2.OutputName = {'y' 'x1' 'x2'};
gff2.InputName = 'd';
gff2.OutputName = 'uffin';
sensor = sumblk('ymeas = y + d');
Ki     = ss(K_lqi2(3),'InputName','ei','OutputName','ui');
Kstate = ss(K_lqi2(1:2),'InputName',{'x1' 'x2'},'Outputname','ustate');
Kfb    = tf(realp('Kfb',1),'InputName','ufbin','OutputName','ufb');
Kff    = tf(realp('Kff',1),'InputName','uffin','OutputName','uff');
sys = connect(...
    sumblk('e = r - ymeas'), ...
    tf(1,[1 0],'InputName','e','OutputName','ei'), ...
    Ki, ...
    Kstate, ...
    plant2, ...
    sensor, ...
    gff2, ...
    sumblk('ufbin = -ui - ustate'), ...    % Note polarity on ui !!!
    sumblk('u     = uff + ufb'), ...
    Kfb, ...
    Kff, ...
    {'r' 'd'}, ...
    {'y' 'ymeas' 'u' 'uff' 'ufb'});

[ymeasconnect,t]   = step(sampleBlock(sys({'ymeas'},'d'),'Kfb',1,'Kff',1));

% H2(s) = ( F(s)*P2(s) + 1*(1 + Kvel*P2(s)*s + Ktheta*P2(s) ) / (1 + Ki*P2(s)/s + Kvel*P2(s)*s + Ktheta*P2(s) )
F      = gff2;
P2     = plant2('y','u');
Ktheta =  K_lqi2(1);
Kvel   =  K_lqi2(2);
Ki     = -K_lqi2(3);  % negative because block diagram has u = ui - ufb, but connect has u = -ui -ufb
s      = tf('s');
H2     = ( F*P2 + 1*(1 + Kvel*P2*s + Ktheta*P2 )) / (1 + Ki*P2/s + Kvel*P2*s + Ktheta*P2);

ymeasalgebra = step(H2,t);

figure
plot(t,ymeasconnect,t,ymeasalgebra)

[John]

I see -- I had missed that you used a different block diagram junction that i used, I think due to how LQI() computes Ki gain (ie returns "-" gain, as you uncovered in your post here):

At the junction between Ki*xi and u_pre, I have a "-,-" for the Ki and sum(Ki*xi) junction block, vs "+,-" in yours, perhaps for the LQI() gain reason.

That would change the back path gains of the outer loop, then (ie the outer loop has an extra -1x)

Originally when building H(s), i had

loop path gains: -Ki*P/s, -Kth*P, -Kw*P*s,

characteristic eqtn: 1 - (Ki*P/s + P*Kth + P*Kw*s),

cofactor / independent loop paths: 1 - (-Kth*P -Kw*P*s), 1-(0)

which led to the same H(s) that you showed.

The difference in junction sign would mean the loop path gains change to +Ki*P/s, -Kth*P, -Kw*P*s, vs (ie change to +,-,-, from -,-,-), so the characteristic eqtn would also change.

If so, that seems to lead to an oddly behaving H(s) (eg the starting value is 0, vs 1), so that doesn't seem to be the right change. Would your H change in the way I described, if the returned Ki is negative from LQI()?

Note: Perhaps i'm still misunderstanding the Ki sign question -- I can confirm that my loop below behaves as expected with that Ki junction as it is now, when assuming that LQI() returns (+,+,-) for (Kth, Kw, Ki). The behavior of the sim, and all Matlab lsims() when using the correspnding state equations, from all seems to behave as expected -- just H(s) does not, however i check it.

[Paul]

Refering to block diagram immediately above:

All I did was change sign of Ki1 AND the sign on the sum junction for Ki*xi. So there should be no change in the overall system.

If P(s) = Y(s)/U(s), then the characteristic equation is:

Del(s) = 1 - Ki1*P(s)/s + Kw*s*P(s) + Kth*P(s)

The forward paths are: Gff2(s) + 1*(1 + Kw*s*P(s) + Kth*P(s))

The closed loop transfer function is:

H(s) = Yafterdis(s) / D(s)

= ( Gff2(s) + 1 + Kw*s*P(s) + Kth*P(s) ) / ( 1 - Ki1*P(s)/s + Kw*s*P(s) + Kth*P(s) )

If I change the sign on Ki1 AND the sign at that sum junction, then the H(s) would be

H(s) = ( Gff2(s) + 1 + Kw*s*P(s) + Kth*P(s) ) / ( 1 + (-Ki1)*P(s)/s + Kw*s*P(s) + Kth*P(s) )

which is mathematically equivalent.

[John]

Thanks. The new H i've been using matches what you just obtained.

I see an odd output when I'd expected a flat line (see below).

I'm using these outputs of H(s), automatically created symbolically based on the loop structure:

H_characteristic = Kth*P + Kw*P*s - (Ki*P)/s + 1 %This matches

H_symbolic = (F*P + Kth*P + Kw*P*s + 1)/(Kth*P + Kw*P*s - (Ki*P)/s + 1) % This also matches

Substituting F by defining first finding a symbolic F using solve(Heqn, F), where

Heqn = H == 0

gives:

F = - Kth - Kw*s - 1/P % Still the same (as expected)

Then use

H_symbolic = subs(H_symbolic)

to update the expression, which gives:

H_symbolic = (Kth*P - P*(Kth + Kw*s + 1/P) + Kw*P*s + 1) / (Kth*P + Kw*P*s - (Ki*P)/s + 1)

This is just 0 / (Kth*P + Kw*P*s - (Ki*P)/s + 1) = 0

1) So, H --> 0 is expected, but in particular I'd also expect h(t) = 0 for all t. But there's an odd response:

Is this numerical inaccuracy? Should this not give a flat line at 0, given H?

2) In case it matters, the code I use to generate the plot reflects the above matlab-found symbolic expression:

F = zpk( -1/P - K(1) - K(2) * s);

H = zpk((F*P + K(1)*P + K(2)*P*s + 1) /...

(-(K(3)*P)/s + K(1)*P+ K(2)*P*s + 1));

And this gives me (for this example)

H =

1.1586e-13 s (s-451.6) (s^2 + 240s + 1.421e05)^4

------------------------------------------------

(s+2019) (s^2 + 240s + 1.421e05)^5

Again, why is this not returning 0?

3) In any case, I added a pair of fast poles to F, to see if the reponse is now matching what's expected in a physical implementation (ie wouldn't expect h(t) to be 0 for all t)

H =

(s+4.155e06) (s+1.713e06) (s+2e05) (s+0.001807) (s-0.001805) (s^2 + 240s + 1.421e05)^5

--------------------------------------------------------------------------------------

(s+1e05)^2 (s+1.71e06) (s+4.156e06) (s+2022) (s^2 + 240s + 1.421e05)^5

Then, step(H) yields

And bode(H) yields

Here, I'm surprised at the phase -- this seems to indicate an already unstable dsiturbance response system. Can this bode plot be correct, or is there some impact frmo similar potential numerical issues from (1) and (2)?

Regardless, the step response doest match the SL model output. So the H construction does work as expected with the Ki sign change, when fast poles are included in F(s).

But since I know the SL model response is stable, that also amplifies the question around the phase.

[Paul]

1) The scale on that plot is 1e-17, so it's basically zero. Can't expect all of those numerical manipulations to form H to numerically result in 0.

2) Same as (1). The gain is 1e-13.

3) That's the Bode plot of H(s), which is the closed loop transfer function, which is not normally used to assess stability. Clearly the system is stable based on the poles of H(s). Not sure why you think that Bode plot indicates otherwise.

At low frequency, the gain is headed towards zero and the phase is 180. That implies the step response will approach zero from below, which is exactly what you see in the time response. Assuming that H was typed in correctly, a negative final value of the step response is to be expected because of the (s-0.001805) term in the numerator.

I can't recreate your results, doing as best I can with what you showed above.

load FF_check
% H2(s) = ( F(s)*P2(s) + 1*(1 + Kvel*P2(s)*s + Ktheta*P2(s) ) / (1 + Ki*P2(s)/s + Kvel*P2(s)*s + Ktheta*P2(s) )
%F      = gff2;
P2     = plant2('theta','u');
Ktheta =  K_lqi2(1);
Kvel   =  K_lqi2(2);
Ki     =  K_lqi2(3);  % use LQI configuration
s      = tf('s');
F      = zpk(-1/P2 - Ktheta - Kvel*s);
H2     = zpk( (F*P2 + Ktheta*P2 + Kvel*P2*s + 1) / (-Ki*P2/s + Ktheta*P2 + Kvel*P2*s + 1) )

Warning: The frequency response has poor relative accuracy. This may be because the response is nearly zero or infinite at all frequencies, or because the state-space realization is ill conditioned. Use the "prescale" command to investigate further.

H2 =
 
  -1.1532e-12 s (s-52.3) (s^2 + 240s + 1.421e05)^4
  ------------------------------------------------
         (s+2019) (s^2 + 240s + 1.421e05)^5
 
Continuous-time zero/pole/gain model.

%{
% Your result
H =
   1.1586e-13 s (s-451.6) (s^2 + 240s + 1.421e05)^4
  ------------------------------------------------
         (s+2019) (s^2 + 240s + 1.421e05)^5
%}

figure
[y1,t1] = step(H2);
plot(t1,y1);

If I were to try to come up with H2 algebraically, I would try something like this:

P2     = zpk(plant2('theta','u'));
Ktheta =  K_lqi2(1);
Kvel   =  K_lqi2(2);
Ki     =  K_lqi2(3);  % use LQI configuration
s      = zpk('s');
F      = (-1/P2 - Ktheta - Kvel*s)

F =
 
  From input "theta" to output:
  -7.0362e-06 (s+4.155e06) (s+1.713e06)
 
Continuous-time zero/pole/gain model.

Define two terms that sum together to form the numerator of H2

numterm2    = 1;
numterm1    = P2*(F + Ktheta + Kvel*s)

numterm1 =
 
  From input "theta" to output "theta":
  - (s^2 + 240s + 1.421e05)
  -------------------------
   (s^2 + 240s + 1.421e05)
 
Continuous-time zero/pole/gain model.

[z,p,k] = zpkdata(numterm1);
format long e
z{1},p{1}

ans = 
     -1.200000000029235e+02 + 3.573825772130968e+02i
     -1.200000000029235e+02 - 3.573825772130968e+02i

ans = 
     -1.200000000000000e+02 + 3.573825728483228e+02i
     -1.200000000000000e+02 - 3.573825728483228e+02i

Obviously (to me) those poles and zeros are supposed to cancel.

Try pole-zero cancellation with default tolerance

numterm1 = minreal(numterm1)

numterm1 =
 
  From input "theta" to output "theta":
  - (s^2 + 240s + 1.421e05)
  -------------------------
   (s^2 + 240s + 1.421e05)
 
Continuous-time zero/pole/gain model.

No change. I recall from some other question that I answered that minreal uses a very tight tolerance for poles and zeros that are close to the imaginary axis. Use a large tolerance to force the cancellation.

numterm1 = minreal(numterm1,1e-6)

numterm1 =
 
  From input "theta" to output "theta":
  -1
 
Static gain.

but numterm1 isn't exactly -1

[z,p,k] = zpkdata(numterm1);
k

k = 
    -1.000000000003124e+00

num    = numterm1 + numterm2;

I'd form the denominator this way:

den = P2*(-Ki/s + Ktheta + Kvel*s) + 1;

Then H2 is

H2     = (numterm1 + numterm2) / den

H2 =
 
  From input "theta" to output "theta":
  -3.1242e-12 s (s^2 + 240s + 1.421e05)
  -------------------------------------
    (s+4.156e06) (s+1.71e06) (s+2022)
 
Continuous-time zero/pole/gain model.

[y2,t2] = step(H2);

H2 is not exactly zero because numterm1 is not exactly -1.

We can make it exactly zero by forcing the value of numterm1

(-1 + numterm2) / (-Ki*P2/s + Ktheta*P2 + Kvel*P2*s + 1)

ans =
 
  0
 
Static gain.

Try with connect()

plant2.C = [plant2.C;eye(2)];
plant2.OutputName = {'y' 'x1' 'x2'};

gff2 = F;
gff2.InputName = 'd';
gff2.OutputName = 'uffin';
sensor = sumblk('ymeas = y + d');
Ki     = ss(K_lqi2(3),'InputName','ei','OutputName','ui');
Kstate = ss(K_lqi2(1:2),'InputName',{'x1' 'x2'},'Outputname','ustate');
Kfb    = tf(1,1,'InputName','ufbin','OutputName','ufb');
Kff    = tf(1,1,'InputName','uffin','OutputName','uff');
sys = connect(...
    sumblk('e = r - ymeas'), ...
    tf(1,[1 0],'InputName','e','OutputName','ei'), ...
    Ki, ...
    Kstate, ...
    plant2, ...
    sensor, ...
    gff2, ...
    sumblk('ufbin = -ui - ustate'), ...    % Note polarity on ui !!!
    sumblk('u     = uff + ufb'), ...
    Kfb, ...
    Kff, ...
    {'d'}, ...
    {'ymeas'});

Warning: The following block inputs are not used: r.

sys = minreal(sys);

3 states removed.

zpk(sys)

ans =
 
  From input "d" to output "ymeas":
  -1.3323e-15 s (s+3.309e05) (s+5.157e06)
  ---------------------------------------
     (s+4.156e06) (s+1.71e06) (s+2022)
 
Continuous-time zero/pole/gain model.

[ymeas,t]   = step(sys({'ymeas'},'d'));
figure
subplot(311),plot(t,ymeas),grid
subplot(312),plot(t1,y1),grid
subplot(313),plot(t2,y2),grid

All three cases are effectively zero and illustrate the numerical complications of evaluating block diagrams that should mathematically result in a zero response.

[John]

"The scale on that plot is 1e-17, so it's basically zero. Can't expect all of those numerical manipulations to form H to numerically result in 0. ... Same as (1). The gain is 1e-13."

I see. I would have thought that the 1e-13 is relatively large for 64-bit operations (precision of 16 decimal digits, so 1000x smaller than this). But, understood that this seems expected :)

"That's the Bode plot of H(s), which is the closed loop transfer function, which is not normally used to assess stability. Clearly the system is stable based on the poles of H(s). Not sure why you think that Bode plot indicates otherwise."

It seems to me that it'd represent the physical response, ie the system's response to a disturbance. Near DC, i'd expect the highest attenuation (smallest impact of D(s)), and consider F somewhat effective up to the crossover frequency. Similarly, F should not signifcantly amplify noise (eg sensors) above crossover, and ideally attenuate as well -- so that brings up that perhaps F would benefit from rolloff at higher freqencies, which is not how it currently is. (To create that rolloff, it's not just a LPF given the FB loop, but haven't looked yet into how to roll off the bode desponse).

For the phase portion: i would have thought H(s) phase margin indicates how robust the design is, eg a small phase margin would indicate poor tolerance to additional delays, modeling imperfections, etc. So at DC, I would expect a large phase margin, but H(s) currently has 180 (= -180, ie ~0 phase margin)

"At low frequency, the gain is headed towards zero and the phase is 180. That implies the step response will approach zero from below, which is exactly what you see in the time response."

Interesting -- why does the phase indicate direction of approach? It should just indicate the phase of the response but not the path, no? Per above, i would have thought that indicates robustness of the loop.

"Assuming that H was typed in correctly, a negative final value of the step response is to be expected because of the (s-0.001805) term in the numerator"

dcgain(H) = -3.2322e-14, which should be considered 0 per the discussion above on resolution. I used the H with fast poles (from my last post):

(s+4.155e06) (s+1.713e06) (s+2e05) (s+0.001807) (s-0.001805) (s^2 + 240s + 1.421e05)^5

--------------------------------------------------------------------------------------

(s+1e05)^2 (s+1.71e06) (s+4.156e06) (s+2022) (s^2 + 240s + 1.421e05)^5

To check, you're looking specifically at (s-0.001805) because it's the only RHP zero; so you're saying there will be an undershoot and slower response because of it. Is this right, and if so, wouldn't the small-magnitude result of dcgain(H) indicate 0, vs a negative value?

"Try pole-zero cancellation with default tolerance

<...>

No change. I recall from some other question that I answered that minreal uses a very tight tolerance for poles and zeros that are close to the imaginary axis. Use a large tolerance to force the cancellation."

Ah -- right; good point and reminder about the tolerance for cancellation. Thanks, that makes sense :)

[Paul]

My only claim about that right-half plane zero is that the final value of the step response would be a negative value, small as it may be for the values of the coefficents as they turned out after all fo the numerical computations.

I'm going to retract my statement about the low frequency phase of H(s) indicating the response will approach it's final value from below. Not sure what I was thinking there.

As for "i would have thought H(s) phase margin indicates how robust the design is....", again, H(s) is the closed-loop system transfer function, which is not used to assess stability margins. For that, you'd need a transfer function of an open-loop system computed at a loop break point for where the stability margins are to be computed.

[John]

"My only claim about that right-half plane zero is that the final value of the step response would be a negative value, small as it may be for the values of the coefficents as they turned out after all fo the numerical computations."

Got it -- thanks. Since it was the same mangitude as the signals considered as 0, I hadn't been assigning meaning to the sign of the ~1e-13 value. But that makes sense to do so.

"H(s) is the closed-loop system transfer function, which is not used to assess stability margins."

True... Not sure why I thought the phase has stability meaning for H; agreed with what you're saying. There should still be some physical significance to the phase (since there's a signal propagating through the closed loop), but i'll think more on the implication separately.

As a further test of F(s) rejection of D(s), I implemented FB+FF in discrete time using Matlab and Simulink:

• discrete gains via discrete lqi() form using the loop time
• discretized F(s) using the loop time with c2d()
• still using a continuous plant (ie digital control of a physical plant).
• Simulink using rate transition blocks for the states x and output y, and a discrete-time integrator. (So, x1 (theta) and x2 (vel) are still continuous-time)
• The discrete-time loop time is fast enough that doubling the rate has no impact on Ym response.

FB and FF control signals both track between the continuous and discrete versions of the system; with the FB loop closed (ie FF not active), the Ym responses match.

But, when adding FF in and checking Ym from FF+FB, the responses do not match: the discrete response substantially differs. This is surprising, as I would have thought ~equivalent command signals would yield ~equivalant responses. Is there any issue with the loop structure below? Perhaps a 1/Ts gain or similar is needed somewhere, eg the integrator, but if so, not sure why.

Here's the system (not attaching a model, just a screenshot). "d" denotes the discrete version.

Here's the continuous version (from the thread above) just for easier comparison; pardon the variable names not matching, but the structures match.

Some plots of this: for a unit step, u_FB and u_FF match for both cts and disc cases. In the plot below only FB is connected to the loop, but i'm showing u_FF as well, to show u_FF seem to match even if they're not applied to the plant.

Zoom:

Then i connect u_FF (FF active). The discrete version has an unexpdected excursion to -400. I'm showing Ym and Ym_d in two plots, for easier viewing due to their scale difference).

Per the cursor, u_FF and u_FF_d version seemingly correctly at ~ -5e7, per the thread above. Also Ym and Ym_d both settle to 0 (Ym_d happens to be at -1 at the cursor).

Even during this unexpected large excursion, u_FB and u_FF both seem to be ~identical (as expected) between the cts and disc versions.

One last test, with FF-only (FB not active). Here, Ym responses match exactly again. This test uses the F(s) meant for FB (-1/P - Kth - Kw*s), ie since FB is intentionally inactive, the response is subpar as expected.

How can it be that all signals match, and responses are the same by themselves, but there is drastically different response with FF+FB?

Simulink solver is set to fixed-step 1e-8, which seems like it should be fast enough.

There doesn't seem to be anything unstable with the discrete version of F(s):

Fd =

-70362 (z+207.6) (z+0.3636)

---------------------------

(z-0.3679)^2

Sample time: 1e-05 seconds

Discrete-time zero/pole/gain model.

And bode of both F(s) and Fd(s) also match:

[Paul]

I might not be following all of the above, but ....

I used Fd as defined above, with a Sample Time of 1e-5. I used the K_lqi2 as have been defined throughout this thread, I replaced the Integrator with a Discrete Time Intergrator, and used ZOH blocks on the feedback paths with a Sample time of 1e-8. This hybrid system yielded the same result for ymeas as the fully continous system, with both FF and FB in use.

From what I gather, you redesigned the LQI gains for the discrete control; are you sure discrete-time controller designed this way would yield a similar response to the continous system. Also, Fd has a Sample Time of 1e-5 as do all of your feedback loops. Ts = 1e-5 would be too large to replicate the dynamics of the continuous-time system, which is why I used 1e-8.

Not sure if this helps ....

[John]

@Paul thanks.

"This hybrid system yielded the same result for ymeas as the fully continous system, with both FF and FB in use."

Huh, interesting that it worked for you. I'll investigate more why the difference in sampling rates matters, as (more below) i thought they should match given c2d() converting an entire system's individual components.

When i swap in Fd(s) (the c2d() version of F(s)) into the continuous model, that model also works for FB+FF. So, Fd(s) seems like it works, so not sure why it doesn't work with the discrete FB controller -- especially when that discrete FB controller also works by itself:

"redesigned the LQI gains for the discrete contro... are you sure discrete-time controller designed this way would yield a similar response to the continous system."

Yes; i've used that successfully several times before. It seems to work in this instance too: the first plot above shows the FB response, with both cts and disc gains. Both match for essentially any loop rate, so it's working. Similarly, FF only also shows the exact same response between disc and cts (per plots above).

But, it's a diffnon-matching response with FF+FB: the disc system has large excursion to -400.

"Also, Fd has a Sample Time of 1e-5 as do all of your feedback loops. Ts = 1e-5 would be too large to replicate the dynamics of the continuous-time system, which is why I used 1e-8."

Hmm, but if this were the cause, wouldn't the FF-only responses also differ? They match between disc and cts.

The loop is running at a given rate (1e-5), so I'd sampled all blocks for that given rate. Is this the wrong utilization of c2d()?

[John]

On further examination of the issue, I think i have a better understanding of what's causing this: assumptions on H(s) that don't seem to hold in the discrete case where there's feedback loop delay, and uFF's impact isn't fed back fast enough via FB's Kth+Kw*s. I don't know what the solution is, though.

The issue seems connected to my question on uFB "mirroring/canceling" uFF, ie lots of FF and FB effort that cancel out to a small number. @Paul you had a good derivation in this post about how this is expected, since the input to plant equation in steady state (in response to disturbance d) is

u = GFF2*d - Kvel*thetadot - Ktheta*theta + Ki*ei --> -1 = Gff + (Kth + Ki*ei), so

"... in steady state

-1 = FF + FB

In the absence of feedforward control, the steady state feedback would be FB = -1

With the feedforward control, we'd have FB = -1 - FF

which is essentially FB = -FF as you've observed because FF is enormous compared to unity for the choice of Gff2."

In the case of a continuous (ie fast feedback) system, this FB ~= - FF holds, and so it's okay if both FF and FB are very large, since they mirror each other and essentially cancel out. H(s) still holds since uFB quickly compensates for the very large uFF. So gff2 = -1/P - Kth - Kw*s has a direct and instantaneous counterpart of +Kth+Kw*s coming from FB, given how H(s) includes both the FB and FF terms.

Essentially: gff2 = -1/P - Kth - Kw*s is valid only when the rest of H(s) is valid since that's how it was derived, ie only valid if FB includes + Kth + Kw*s.

But in the discrete case, FB doesn't immediately include FF's impact via + Kth + K*s , since a loop iteration time is required before the impact shows up in FB, ie one loop delay before FB can mirror/cancel the majority of uFF's impact.

Even at high discrete-time loop rates relative to the plant time constant (~1000x), given the very large FF command driving the plant, the plant changes state 100s or 1000s times larger than the input disturbance (eg -400 in the case above).

As the loop rate increases to 10^5 - 10^6x the plant time constant -- eg your above test of 10^8 Hz loop rate -- the loop essentially becomes continuous, so once again H(s) --> 0 with gff2(s).

This time delay between FF and FB seems a fundamental issue with state feedback and disturbance FF on practical discrete-time systems, given how F(s) (gff2) was derived. Of course, the practical input command limit will change the size of the excursion when FF and FB are off by one cycle, but a complete discrete-time base approach should also take into account the discrete case (vs relying on practical limits).

Overall, not sure how to address this. Slowing down gff2 evolution with slower poles won't help: Y is never bounded to [0,1] unless the loop rate is basically infinite (ie cts), and that's impractical, so it's worse to include FF than to not even use it.

Something is needed like solving for F when H includes both step k and k-1, so perhaps converting H(s) to H(z^-1) and seeing what comes out of that. But i think that wouldn't help...