This note describes how to design a PID controller for a system defined by second order differential equation based on requirements for a step response specified by the rise time and the settling time.
The goal is to determine the three PID parameters (\(K_{i},K_{d},K_{p}\)) from the plant transfer function and \(t_{r},t_{s}\) (rise time and settling time).
Consider the following mechanical system
\(M\) is the mass of the car, \(b\) is the damping coefficient and \(k\) is the spring constant. To illustrate, assuming standard SI units:\begin {align*} M & =1\ kg\\ b & =10\ Ns/m\\ k & =20\ N/m \end {align*}
The first step is to derive the mathematical model for the system. This means finding a differential equation that relates the output (the displacment \(y(t)\)) to the input, which is the applied force \(u(t)\).The fFriction force between the mass M and the ground is ignored in this example.
The first step is to make a free body diagram
Applying Netwon laws gives \begin {align*} F & =ma\\ u(t)-by^{\prime }(t)-ky(t) & =My^{\prime \prime }(t) \end {align*}
or\[ My^{\prime \prime }(t)+by^{\prime }(t)+ky(t)=u(t) \] Taking Laplace transform and assuming zero initial conditions gives \[ Ms^{2}Y(s)+bsY(s)+kY(s)=U(s) \] The transfer function is defined as the ratio of the output to the input in the Laplace domain. Here the input is \(u(t)\), which is the external force, and the output is \(y(t)\) which is the displacement. Taking the Laplace transform of the above differential equation gives the transfer function \[ sys(s)=\frac {Y(s)}{U(s)}=\frac {1}{Ms^{2}+bs+k}\] Using block diagram the transfer function is illustrated as
The PID controller is now added. The transfer function of the PID controller itself is\[ K_{p} + K_{d} s + \frac {K_{i}}{s} \] The controller is added to the system and the loop is closed. The following diagram represents the updated system with the controller in place
Let \(L(s)\) be the open loop transfer function \[ L(s) = \frac {K_{p}+K_{d} s + \frac {K_{i}}{s}} {M s^{2} + b s + k} = \frac {K_{p} s + K_{d} s^{2} + K_{i}} {M s^{3} + b s^{2} + k s} \] Hence the closed loop transfer function is\[ G\left ( s\right ) =\frac {L\left ( s\right ) }{1+L\left ( s\right ) }=\frac {K_{p}s+K_{d}s^{2}+K_{i}}{Ms^{3}+bs^{2}+ks+K_{p}s+K_{d}s^{2}+K_{i}}\] Therefore \begin {equation} G\left ( s\right ) =\frac {\frac {1}{M}\left ( K_{d}s^{2}+K_{p}s+K_{i}\right ) }{s^{3}+s^{2}\left ( \frac {b+K_{d}}{M}\right ) +s\left ( \frac {k+K_{p}}{M}\right ) +\frac {K_{i}}{M}} \tag {1} \end {equation}
The closed loop transfer function (1) shows there are three poles.
Putting one pole at a distance of \(5\zeta \omega _{n}\) away from the imaginary axis, while the remaining two poles are the dominant poles results in the following diagram
The denominator of equation (1) can be rewritten as\begin {align*} s^{3}+s^{2}\left ( \frac {b+K_{d}}{M}\right ) +s\left ( \frac {k+K_{p}}{M}\right ) +\frac {K_{i}}{M} & \equiv \left ( s+5\zeta \omega _{n}\right ) \left ( s^{2}+2\zeta \omega _{n}s+\omega _{n}^{2}\right ) \\ & \equiv s^{3}+2\zeta \omega _{n}s^{2}+s\omega _{n}^{2}+5\zeta \omega _{n}s^{2}+10\zeta ^{2}\omega _{n}^{2}s+5\zeta \omega _{n}^{3}\\ & \equiv s^{3}+s^{2}\left ( 7\zeta \omega _{n}\right ) +s\left ( 10\zeta ^{2}\omega _{n}^{2}+\omega _{n}^{2}\right ) +5\zeta \omega _{n}^{3} \end {align*}
Equating coefficients gives \begin {align*} \frac {b+K_{d}}{M} & =7\zeta \omega _{n}\\ \frac {k+K_{p}}{M} & =10\zeta ^{2}\omega _{n}^{2}+\omega _{n}^{2}\\ \frac {K_{i}}{M} & =5\zeta \omega _{n}^{3} \end {align*}
Solving for PID parameters results in \begin {align} K_{d} & =7M\zeta \omega _{n}-b\tag {2}\\ K_{p} & =M\left ( 10\zeta ^{2}\omega _{n}^{2}+\omega _{n}^{2}\right ) -k\nonumber \\ K_{i} & =5M\zeta \omega _{n}^{3}\nonumber \end {align}
These are the PID parameters as a function of \(\omega _{n}\) and \(\zeta \).
\(\omega _{n}\) and \(\zeta \) are determined in order to obtain the PID parameters.
The time response specifications are now introduced in order to determine these parameters. Assuming these are the time domain requirments
Using the following for \(2\%\) criterion
\begin {equation} t_{s} = \frac {4}{\zeta \omega _{n}} \tag {3} \end {equation} And the rise time is given by\[ t_{r} = \frac {\pi -\beta }{\omega _{d}} \]
But \(\beta =\arctan \left ( \frac {\omega _{n}\sqrt {1-\zeta ^{2}}}{\zeta \omega _{n}}\right ) \), hence \begin {equation} t_{r}=\frac {\pi -\arctan \left ( \frac {\sqrt {1-\zeta ^{2}}}{\zeta }\right ) }{\omega _{n}\sqrt {1-\zeta ^{2}}} \tag {4} \end {equation} From (3) and (4) \(\omega _{n},\zeta \) are solved for \begin {align*} 2 & =\frac {4}{\zeta \omega _{n}}\\ 0.1 & =\frac {\pi -\arctan \left ( \frac {\sqrt {1-\zeta ^{2}}}{\zeta }\right ) }{\omega _{n}\sqrt {1-\zeta ^{2}}} \end {align*}
Solving numerically gives
Hence the solution is \[ \zeta =0.1176 \] and \[ \omega _{n}=17\ \ \text {rad/sec}\] Substituting the values for \(\omega _{n}\) and \(\zeta \) in (2), and the values given for \(M,b\) and \(k\), gives the PID parameters \begin {align*} K_{d} & =7(1)(0.1176)(17)-10\\ & =3.9944\\ K_{p} & =1(10(0.1176^{2})(17^{2})+(17)^{2})-20\\ & =308.97\\ K_{i} & =5(1)(0.1176)(17)^{3}\\ & =2888.8 \end {align*}
Using Matlab, the step response is found