M-file PID control of mass-spring-damper

[Barry Dolan]

Hi,

I am trying to implement a PID controller on a mass-spring-damper system in continuous time. I would like to do this in a Matlab m-file, rather than in Simulink.

I have checked what the step response of the system looks like in the s-domain (i.e. with the system transfer function together with the built-in Matlab command "step") using the same PID tuning, and I know this answer to be correct.

However, when I simulate the system in the time domain, the response is very different.

I have seen some other examples of PID implemented in an m-file, but these all seem to use a transfer function based approach, and I can't figure out how to get the correct answer using my approach.

I have copy/pasted two m-files below. The first of these defines the model parameters and calls the second m-file function using ode45.

I would greatly appreciate if someone could advise me on this topic.

Thank you.

% Here is the first file:
% This program simulates the step response of a mass-spring-damper system,
% using PID control.

clear all, close all, clc

global M b k Kp Ki Kd x_ref

% Model Parameters
M = 1; b = 10; k = 20;

% The reference position is 1 metre
x_ref = 1;

% PID
Kp = 350; Ki = 300; Kd = 50;

sim_time = 2;

% Initial states
x = 0; dx = 0;
p0 = [x dx]';

options = odeset('RelTol',1e-6);
rhsfun = @(t,p) sos(t,p);
[t,p] = ode45(rhsfun,[0 sim_time],p0,options);

figure('Name','System Response to Step Input','NumberTitle','off')
plot(t,p(:,1)), grid on, xlabel('Time[s]'), ylabel('Amplitude [m]')

% Here is the second file:
function dp = sos(t,p)

global M b k Kp Kd Ki x_ref

persistent ep_old ep ed ei t_old

% Assigning values to persistent variables on first call of function
if isempty(ep_old)
ep_old = x_ref - p(1);
ep = x_ref - p(1);
ed = 0;
ei = 0;
t_old = 0;
end

dt = t - t_old;

% Calculating the PID controlled system input
F = Kp*ep + Kd*ed + Ki*ei;

% State update
dp = [p(2); M\(F - b*p(2) - k*p(1))];

% Calculation of the error and its derivative and integral
% Note: the error is updated only when the solver steps forwards in time
if dt > 0
ep = x_ref - p(1); %Calculating P
ed = (ep - ep_old)/dt; %Calculating D
ei = ei + ep*dt; %Calculating I
ep_old = ep;
end

t_old = t;

[Torsten]

Solving the system
dp1/dt = p2
dp2/dt = (F - b*p2 - k*p1) / M = (Kp*(x_ref-p1) - Kd*p2 + Ki*p3 - b*p2
- k*p1)/M
dp3/dt = x_ref - p1
with p1(0)=p2(0)=p3(0)=0
should work better.
Note that the adpative time stepping of ODE45
does not allow calculation of int_ep and d/dt(e_p)
during the integration as you do above.
You need a solver with fixed time steppin for this
(as supplied in Simulink, I think).

Best wishes
Torsten.

[Barry Dolan]

"Barry Dolan" wrote in message <im9vqs$8jt$1...@fred.mathworks.com>...

Thanks for the advice Torsten, I have implemented your suggestion and included the error as a model state.

However still I observe a discrepancy between the response I calculate in the frequency domain using PID control compared with the time domain response using PID control.

If I compare the two solutions using any of P, PI or PD control, there is no discrepancy between the frequency and time domain step responses.

I wonder if there is some obvious reason for this which I am overlooking?

Bests,

Barry.

[Torsten]

On 22 Mrz., 13:50, "Barry Dolan" <barry...@gmail.com> wrote:
> "Barry Dolan" wrote in message <im9vqs$8j...@fred.mathworks.com>...

> Barry.- Zitierten Text ausblenden -
>
> - Zitierten Text anzeigen -

Sorry, but I'm not used to your vocabulary.
How do you get the response in the frequency domain and compare it
with the time domain ?

Best wishes
Torsten.

[Barry Dolan]

Torsten <Torsten...@umsicht.fraunhofer.de> wrote in message <e5979295-09b7-4593...@s9g2000yqm.googlegroups.com>...

Yes, my vocabulary was incorrect, sorry. I should have included the m-file to illustrate my problem. Below is the code I am using to generate the step response using the transfer function (which I incorrectly said was a "frequency domain response") and Matlab command "step". The step response arising from this method is not in agreement with that arising from the code which I have posted at the beginning of this thread (using ode45). I hope that now my problem is clearer to you.

% PID Controller
% Closed-loop TF: X(s)/F(s) = (Kd*s^2 + Kp*s + Ki)/(s^3 + (10 + Kd)*s^2 + (20 + % Kp)*s + Ki)
Kp = 350; Ki = 300; Kd = 50; num = [Kd Kp Ki]; den = [1 10+Kd 20+Kp Ki];
t = 0:.01:2;
figure('Name','System Response to Step Input with PID Controller','NumberTitle','off')
step(num,den,t)