Discretized transfer function feedback implementation

[jblc]

I'm applying feedback to a discretized time-series transfer function in software, to simulate how a physical system moves (there are other ways, I'm doing it this way for a reason).
The step response of the derived time-series system IS what i'd expect when I match it with a matlab "step" command's output (applied to the original laplace transfer function), but is NOT what I'd expect when I try closing a loop and adding feedback.
Is there something fundamentally flawed with the proportional feedback approach as shown below?

Say I have a continuous open-loop laplace-domain transfer function of a system, 'TF(s)'. Discretize this to z-domain, 'TF(z)'.
Create a time-series from this, with the output 'x' as a function of the input 'u'. I won't show the actual numbers, but the form is x(k+1) = ... *x(k) + ...*x(k-1) ... *u(k) + ...*u(k-1)... etc

The unit-step response is this:

FOR k = 1 to whatever time step
 u(k+1) = 1 <--unit step input command
 x(k+1) = ... *x(k) + ...*u(k) ... etc
 housekeeping: shift buffer of x and u vectors, etc.
END

The above appears to work, so I try a feedback loop to track a reference input 'ref':

ref = 1;
K = 1; <- feedback constant
FOR k = 1 to whatever time step
 Error = ref - x(k)
[b] u(k)[/b] = K * Error <- input depends on error; basic proportional feedback. [u]Use this u(k) in next line.[/u]
 x(k+1) = ... *x(k) + ...*[b]u(k)[/b] ... etc
 housekeeping: shift buffer of x and u vectors, etc.
END

The latter code doesn't oscillate/settle around 'ref' like I'd expect from simple error feedback. As time grows, x(k) settles on a random values depending on the combination of K and 'ref'. Eg K=1,ref=2 settles at 1. K=2,ref = 0.5 settles at 0.666. K=10,ref=0.5 settles at 3.33.
Obviously, Normally, 'x' would settle around 'ref' with its oscillation freq dependent on K.

Thoughts?

 

[PNachtwey]

K is really your proportional gain.
I would need to see the coefficients to be able to determine what is really going on but K is probably too big. It should be easy enough to calculate a proportional gain Kp that will not oscillate if you have the open loop coefficients for the difference equation. x should settle around the ref if the open loop transfer function has a gain of one. From what you write I would say it doesn't. What are the coefficients?




Peter Nachtwey
Delta Computer Systems

http://www.deltamotion.com

 

[jblc]

Thanks for the reply. Yes, that's right, K is the feedback gain.

I don't believe K is too large: the settling value changes as K varies, but the waveform never actually changes form; it always oscillates. Agreed, small enough K should cause the closed-loop system response to be massively underdamped and not reach its final ref value. For some reason, I never see that in the feedback implementation written above.

1. The example function used is standard 2nd order:
TF = wn^2 / [s^2 + 2*zeta*wn*s + wn^2]
with wn = 1884 rad/s (300 Hz), and zeta (damping factor) = 0.4.
(See the P.S. for the actual discrete transfer function and coefficients.)

2. I think it may be an issue with the discrete implementation or my understanding of what's going on, and I'm curious if my methodology in this part of it is suspect somehow, or if it seems sounds. With this TF, as the freq of input u goes towards zero , s-->0. Thus TF-->1 via wn^2/wn^2=1. So that is consistent with what I see.

Adding feedback should make the settling value --> ref if K is high enough, since CL = K*TF/(1+K*TF). As K gets larger, CL-->1 and tracks ref. I don't know why this isn't happening with the second snippet of pseudo-code in the original post.

P.S. Discrete derivation
The discrete equivalent of TF above is [0.01706*z + 0.01643] / [ z^2 - 1.86*z + 0.8931].
Call 'num' the numerator coefficients = [0, 0.01706, 0.01643], and 'den' the denominator coefficients = [1, - 1.86, 0.8931].
uVect(k=0) = [0 0 0], yVect(k=0) = [0 0 0].
The time series is given by standard form as y(k+1) = [ sum{ num.*uVect(k) } - sum{ den[2:3].*yVect[2:3](k) } ] / den[1].
Then each time step I shift uVect and yVect by one element, so 3 elements of past history are needed in memory.
( A subvector of elements 2 to 3 is selected as [2:3] )

 

[PNachtwey]

That looks like the transfer function for a mass on a spring.
No problem. The document below was done for the plcs.net forum years ago. I wanted to impress upon the forum the advantages of using a derivative gain but this pdf also applies to this problem.
The way I would approach the problem is in pages 1-12. Pandiani, another forum member, took the challenge of trying with just a PI controller only using MatLab in pages 13-14. Pages 15+ are my attempt to control a mass on a spring with a PI controller. I did better than Pandiani but not as good as my solution on pages 1-12.

http://deltamotion.com/peter/Mathcad/Mathcad - T0C1 MassOnASpring-PID.pdf

A P or PI controller will not work well with the non-integrating under damped system. 2 gains are required to place the two complex poles. The integrator in the PID comes with its own poles to altogether there are 3 closed loop poles and 3 gains. On page 12 you can see I can place the 3 poles on the negative real axis in the so domain so there will be NO overshoot.






Peter Nachtwey
Delta Computer Systems

http://www.deltamotion.com

 

[jblc]

Thanks. (That is a very nice and clearly written document, by the way!). Yes, this is a simple spring-damper. I understand in-depth PID (and more sophisticated controllers).

I recognize that two gains are needed for two complex poles given the oscillatory and damped nature of the system, but my explicit purpose for this simulation is to only have a P gain. I realize the output will oscillate, but that's fine: it will oscillate around the setpoint if K is high enough (which it is).

The question is, why is this not oscillating around ref. Matlab's built-in functions work fine as expected with just a K gain, tracking ref. Eg using Matlab's CL = feedback(ref*K*TF, 1) and step(CL) functions will give an appropriate response that is as-expected, as K decreases and grows.
However, the discrete-time version from my original post does not, and it's not clear where the difference comes from...the discrete version should be roughly equivalent to the matlab response.

 

[xnuke]

Please provide a sample period for your discretized TF so I can duplicate it. I'm guessing it's somewhere around 0.0001 based on a few trial and error iterations. I was only able to approximate the coefficients you gave in your post when using the c2d command and I'd like to match your result.

Have you looked at the root locus or Bode plot for your discrete TF? Use MATLAB's sisotool command if you have the Control System Toolbox. Based on my approximation of your discrete TF, it doesn't take much gain (K) before you exceed the gain margin and make your loop unstable. That could be why you're seeing a difference between the continuous time and discrete time simulations. Since your continuous time TF is second-order, it has infinite gain margin. Your discrete time TF does not - gain margin is less than 20 dB in my approximation.

xnuke
"Live and act within the limit of your knowledge and keep expanding it to the limit of your life." Ayn Rand, Atlas Shrugged.
Please see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[jblc]

Yep, exactly on the time. I tested various sample periods, from 1e-4 to 1e-7. The results are almost the same, so I stuck with 1e-4 sec for the numbers given in the post.
I've checked the bode plots and they look fine, when overlaid on top of the continuous bode.

You're right, good call. At 1e-4 the gain margin isn't as high for the discrete system. I just tried 1e-5 which has stability to 300 Hz.

However, the output plots from the time-series is the same for 1e-5 compared to 1e-4, so feel free to use 1e-5 in your testing (I'll use 1e-5 as well; won't match the numbers from the post, but we'll both have the same thing in matlab).

EDIT: Here is my matlab code. Actually, the discrete and continuous versions seem like they're matching now. This looks as expected. Thoughts?

%discrete TF test
clear all; close all;
Ts = 1e-6;
Tmax = 0.01;
tVect = [0:Ts:Tmax];
N = size(tVect,2);
wnHz = 300;
wnRadPs = wnHz * 2 * pi;
zeta = 0.3;

% Plot continuous and discrete TFs
s = tf('s');
tf = wnRadPs^2 / (s^2 + 2 * zeta * wnRadPs * s + wnRadPs^2);
opts = bodeoptions('cstprefs'); opts.FreqUnits = 'Hz'; opts.Grid = 'on';
figure;
bodeplot(tf, opts);
hold on;
tf_disc = c2d(tf, Ts)
bodeplot(tf_disc, opts);
legend('continuous matlab', 'discrete matlab');

% Find discrete coefficients:
[d_t_num d_t_den] = tfdata(tf_disc)
denLen = length(d_t_den{1});
% Prep y and u history vectors
ukvect = zeros(size(d_t_num{1}))';
ykvect = zeros(size(d_t_den{1}))';
yk = zeros(N,1);    % Store state history for plotting

%%%%%%%%%%%%%%%%%%
useFeedback = 1;
%%%%%%%%%%%%%%%%%%
ref = 10;   % Reference (optional based on useFeedback)
Kp = 10;   % Gain
u = 0.1;      % OL input (optional based on useFeedback)
for i = 1:N
    %insert new control signal here based on FB:
    yk(i) = (sum( d_t_num{1}.*ukvect' ) - sum( d_t_den{1}(2:denLen).*ykvect(2:denLen)' ) )/ d_t_den{1}(1);
    ukvect = circshift(ukvect, 1);
    ykvect(1) = yk(i);
    ykvect = circshift(ykvect, 1);
    % Optional feedback.
    if useFeedback == 1
        E = (ref - yk(i));
        ukvect(1) = Kp * E;
    else
        ukvect(1) = u;
    end
end
figure;
hold on;
plot(tVect, yk, '--b')
hold on;
% Compare with ideal matlab CL or OL feedback.
if useFeedback == 0
    % OL step always matches well.
    step(u*tf, Tmax, 'r');
else
    % CL never matches for the discrete version.
    step(ref*feedback(Kp*tf,1), Tmax, 'r');
end
legend('discrete time-series', 'continuous matlab');

 

[xnuke]

I executed your code as written above, first with sample period 1e-6, then 1e-4. For the smaller sampling period, the discrete-time system is stable for a controller gain of 10 since the gain margin is 58.6 dB. For the larger sampling period, it is unstable for a controller gain of 10 with a gain margin of only 18.8 dB. That's why the continuous-time and discrete-time system step responses look the same for the smaller sampling period but not the larger - the larger system is unstable for your chosen proportional gain of 10. I still recommend you import both continuous-time and discrete-time systems into the SISOTool GUI for analysis. You can see the root locus plot, Bode plot, step response plot, etc., by using the GUI. You can then easily see how they all behave as you change the controller gain.

I also see a few things I recommend you change in your code:

You called your continuous-time LTI system object "tf" but you shouldn't do that as "tf" is the name of a native MATLAB CST function. You should always check to make sure you're not naming variables or objects the same name as any native function, because once you declare a variable with the same name, then you can't use that function. Change it to anything other than "tf".

I'd use stairs(tVect, yk, 'b') instead of plot(tVect, yk, '--b') since you're plotting a discrete-time signal (there's no real need to change the line type to dashed if you use stairs).

You're Bode plots look the same for the code above with sample period 1e-6, but they really aren't. Your frequency axis for the code above simply doesn't go high enough. If it went above 10[sup]4[/sup] Hz the two Bode plots would diverge. Try specifying the lower and upper frequency limits in your bodeplot command so that it's wide enough to show the Nyquist frequency and you'll see what I mean.

Apart from writing and implementing the difference equations in code that could be used in a microcontroller, I'm not sure why you aren't simply using the built-in MATLAB step command with your discrete-time system for comparison with the continuous-time system. Care to clarify?

xnuke
"Live and act within the limit of your knowledge and keep expanding it to the limit of your life." Ayn Rand, Atlas Shrugged.
Please see FAQ731-376 for tips on how to make the best use of Eng-Tips.

 

[PNachtwey]

I got a reasonable response with a sample time of 0.0001 seconds. The oscillation settles out after about 7 milliseconds.

http://deltamotion.com/peter/Mathcad/Mathcad - jblc.pdf

I didn't try to optimize the proportional gain because it is futile. With a full PID and using Ackermann's method I can get a critically damped response.

http://deltamotion.com/peter/Mathcad/Mathcad - jblc.pdf

The file is a pdf file. The forum only lets me post pictures.

Peter Nachtwey
Delta Computer Systems

http://www.deltamotion.com