2.1. Control software development#
This lecture will look at the information processing part of a mechatronic system in terms of control loops and digital implementation of controllers.
Information processing part of a mechatronic system:
2.1.1. Digital controller development cycle#
Software components and digital controllers are becoming essential in the development of new automated systems. They offer a wide range of functions at low cost and are flexible in terms of use and maintenance. However, development times are potentially long with specific stages and test times to detect coding errors.
To maximise the benefits of the software and minimise the disadvantages of development, engineers have access to tools at the various stages of the V cycle:
MIL or Model In the Loop: the controller principle is tested by simulation.
RCP or Rapid Control Prototyping: a digital implementation of the controller is implemented in a powerful dedicated processor. The code is generated using specific tools. The physical system is real and can take the form of a reduced-scale model or part of the system.
SIL or Sofwtare In the Loop: the controller software code is computer-tested on a simulated physical system.
HIL or Hardware In the Loop: the software code is implemented in the product’s digital controller. This software/controller assembly is tested electronically on a simulated physical system in an HIL digital target.
Software V design cycle:
Question: Using the example of the control of a gyropode (segway), indicate which stage in the development of the digital controller each image below corresponds to.
HIL quiz:
Answer
(d) MIL: Simulation model with controller principle
(b) RCP: Cylinder test bench with DSPACE controller
(a) SIL: Simulation model with software and clock sampling
(c ) HIL: microcontroller electronic board with HIL target
2.1.2. Continuous controller#
2.1.2.1. Nested loops of a position control#
We will illustrate the effects of digitisation on a launcher nozzle TVC controller. The TVC actuator considered here is an EMA (Electro-Mechanical) linear actuator. Its body is anchored on the booster frame while its rod is attached to the nozzle. In the normal mode, the TVC Actuator must perform a closed loop angular position control of the nozzle, following the demand of the flight control computer and rejecting the disturbances generated by the pulsed jet throught the nozzle. The TVC Actuator must also behave as a damper in order to limit the mechanical resonance effects.
TVC actuator: position and tests
Each of the 2 actuators consists of a brushless permanent magnet synchronous motor that drives a roller screw through a gearbox. This way the rotational motion of the motor is transformed into a linear motion. A typical and generic architecture of small loop controller (big loop = launcher control) is represented below. Classical approach for the design of loop control foresees multiple nested loops (cascade):
The inner current or torque loop is generally associated to a much larger bandwidth than the other loops.
A velocity loop: for EMA, the motor angular position is measured using a resolver and derived to obtain a velocity feedback. A force sensor is added to drive the actuator in order to reduce the transient force and to damp oscillations.
A external position loop: technical constraints impose to have a feedback on actuator position rather than nozzle position itself (which is the controlled parameter). LVDT sensor is generally used to measure the length variation of the actuator.
EMA actuator:
Question: Explain in simplified terms how force feedback can be used to introduce damping.
Damping with forcefeedback
The force feedback is achieved by measuring the force generated by the actuator, using an instrumented rod end for example, and by using this measurment as a setpoint on the speed loop.
A proportional link is therefore created between the actuator force and the linear speed, similar to that of a viscous damper. The gain applied on the force feedback enables to set the damping.
2.1.2.2. Speed loop of an EMA#
We will focus here on the speed loop. We will represent the mechanical assembly nozzle / actuator by the diagram a) below with:
the mass/spring/viscous friction assembly for the nozzle
the rack-and-pinion assembly representing the screw/nut system and the EMA motor
Mechanical models:
Exercice: Explain how to represent the previous model a) by the equivalent model in rotation b). Determine the transfer function relating torque to motor speed under these conditions.
Equivalent parameters
The equivalent parameters must make it possible to represent :
the kinetic energy \(\frac{1}{2}M_{nozzle}v^2=\frac{1}{2}J_{eq}\Omega^2\) thus \(J_{eq,nozzle}=M_{nozzle}p^2\) with \(p\) the pitch in [m/rad]
orelastic energy \(\frac{1}{2}K_{nozzle}x^2=\frac{1}{2}K_{eq}\theta^2\) thus \(K_{eq,nozzle}=K_{nozzle}p^2\) with \(p\) the pitch in [m/rad]
stored in the system as a whole.
Finally:
\(J_{eq}=J_{mot}+M_{nozzle}p^2\)
\(K_{eq}=K_{nozzle}p^2\)
\(f_{eq}=f_{nozzle}p^2\)
Transfer function
\(\frac{\Omega(p)}{T(p)}=\frac{p}{K_{eq}+f_{eq}p+J_{eq}p^2}\)
The Bode diagram of the transfer function can be displayed using the python-control library.
try:
import control as ct
except ImportError as e:
# Install Python-Control package
!pip install control
import control as ct
# Parameters
pitch=1.25e-3 # [m/rad] pitch
Jmot=3e-3 # [kg] motor inertia
fn=5475 # [N/(m/s)] nozzle equivalent viscous friction
Mn=770 # [kg] Nozzle equivalent mass
Kn=478000 # [N/m] Nozzle equivalent stiffness
# Rotational equivalent parameters
Jeq=Mn*pitch**2
Keq=Kn*pitch**2
feq=fn*pitch**2
Jtot=Jmot+Jeq
# Transfer function
sys_meca = ct.tf([1, 0], [Jtot, feq, Keq])
sys_iLoop = ct.tf(1, [1/(2*3.14*500)**2, 2*0.7/(2*3.14*500), 1])
sys_open=sys_iLoop*sys_meca
# Bode diagram
mag, phase, omega= ct.bode_plot(sys_open, plot=True, Hz=True, initial_phase=90,margins=True)
Exercice: Propose a corrector with zero static error and a closed bandwidth of 40 Hz. Make sure you have a minimum phase margin of 45° and a minimum gain margin of 10 dB.
Controller definition
The static error requires the addition of an integral term. We can therefore propose a PI corrector if the stability conditions allow.
\(C(p)=K_p(1+\frac{1}{\tau_i p})=\frac{K_p+K_p\tau_i p}{\tau_i p}\)
The bandwidth in the closed loop can be approximated in the open loop by the frequency with unity gain.
# PI controller
Kp=1.05
tau_i=1/(5*2*3.14)
PI=ct.tf([Kp*tau_i, Kp], [tau_i, 0])
# Open loop
PISys=PI*sys_open
# Bode diagram
mag, phase, omega= ct.bode_plot(PISys, plot=True, Hz=True, margins=True, grid=True, initial_phase=0)
print("PI controller:")
print("Kp = %.2f"% Kp)
print("Tau = %.2e"% tau_i)
print("Ki = Kp/Tau_i = %.2e"% (Kp/tau_i))
PI controller:
Kp = 1.05
Tau = 3.18e-02
Ki = Kp/Tau_i = 3.30e+01
2.1.3. Digitizing a controller and the effect of sampling frequencies#
2.1.3.1. Delay in digital controller#
The previous synthesis of the corrector was done in the continuous domain. However, its implementation will most often be digital rather than analogue. Digital control (using a µC microcontroller or a DSP signal processor) offers the following advantages [Ellis, 2012]:
lower noise sensitivity, better accuracy, less performance drift
ease of changing gains or parameters.
However, digital control has one main weakness: the introduction of delay in the control loop. These delays generate phase shifts that can lead to instability or limit the gains and ultimately the bandwidth of the system.
Digital systems operate in discrete time steps. At regular intervals, the controller reads the feedback and command information, executes control algorithms, and outputs control orders. Between intervals, the output is held at a fixed level. This is what gives digital control waveforms their characteristic “stair step” appearance. The figure below shows a sinusoidal waveform at its sampled equivalent. The effect of sampling with a period \(T_s\) can be approximated as a pure delay of \(T_s/2\).
Sampling effects:
Exercice: Characterise the effect of sampling on gain and phase. Perform the numerical applications for a sampling period of 8 ms (125 Hz).
Note: Delays can be approximated by Padé transfer functions with python implementation in the control package. A first-order Padé approximation of a T time delay is \(e^{-Tp}=\frac{1-\frac{T}{2}p}{1+\frac{T}{2}p}\)
delay=1/125/2
#pade=ct.tf(ct.pade(delay, 1))
num,den=ct.pade(delay, 3)
pade3=ct.tf(num,den)
num,den=ct.pade(delay, 1)
pade1=ct.tf(num,den)
mag1, phase1, omega1= ct.bode_plot(pade1, plot=False, omega_limits=[1,2*3.14*100], Hz=False)
mag3, phase3, omega3= ct.bode_plot(pade3, plot=False, omega_limits=[1,2*3.14*100], Hz=False)
import matplotlib.pyplot as plt
import numpy as np
freq1=omega1/2/np.pi
freq3=omega3/2/np.pi
plt.tight_layout()
fig, (ax1, ax2) = plt.subplots(2) # get subplot axes
plt.sca(ax1) # magnitude plot
plt.plot(freq3, [1]*len(omega3),'r--', label='Exact')
plt.plot(freq1, mag1,'b--', label='Pade 1')
plt.plot(freq3, mag3,'g--', label='Pade 3')
plt.ylabel('Gain')
plt.legend()
plt.sca(ax2) # phase plot
plt.plot(freq3,-delay*omega3*180/np.pi,'r--', label='Exact')
plt.plot(freq1, phase1*180/np.pi,'b--', label='Pade 1')
plt.plot(freq3, phase3*180/np.pi,'g--', label='Pade 3')
plt.xlabel('Frequency (Hz)')
plt.ylabel('Phase (°)')
plt.legend()
<matplotlib.legend.Legend at 0x7f80c39878e0>
<Figure size 640x480 with 0 Axes>
Note: Delays can come from a variety of sources:
Sample-and-hold: 0.5T
Calculation delay: 0.1T-0.9T
Velocity estimation: 0.5T (motion systems)
Their effects are added together to produce a total delay. The anti-aliasing filter can also create phase shifts that need to be taken into account. The sampling frequency must be chosen carefully with all these effects.
Exercice: Evaluate the effect of sampling at a period of 8 ms on the previous corrector. Suggest a sampling frequency compatible with the desired performance.
PISysPade=PISys*pade3
mag, phase, omega= ct.bode_plot(PISysPade, plot=True, Hz=True, initial_phase=0, margins=True)
# find unity gain frequency
difference_array = np.absolute(mag-1)
index_wu = difference_array.argmin()
# plot Nichols around unity gain frequency
w=np.linspace(omega[index_wu]/100,omega[index_wu]*5,100)*2*3.14
ct.nichols(PISysPade, omega=w,grid=True)
Sampling frequency selection
Gain margins and phase margins have fallen sharply compared with the continuous configuration. The sampling frequency must therefore be increased. Assuming a accepted phase variation of 5° for the unity gain angular frequency \(\omega_u\):
\(\omega_u \frac{T_{total}}{2} = \frac{5\pi}{180}\)
\(T_{total} = \frac{10\pi}{180\omega_u}\)
If we assume this total delay \(T_{total} = 1.5 T_{s}\), we obtain the minimum sampling frequency \(f_{s}=\frac{1}{T_{s}}\)
w_u=40*2*3.14 # [rad/s] unity gain angular frequency
Ttotal=10*3.14/180/w_u # [s] total delay
fs=1/(Ttotal/1.5) # [Hz] sampling frequency
print("Sampling minimum frequency = %.2f Hz"%fs)
Sampling minimum frequency = 2160.00 Hz
The following figure shows the effects of sampling on stability trough Modelica simulations. It can be seen that (without taking into account delays due to speed calculation or calculation times) a sampling time of 1 ms gives a behaviour identical to that of the continuous controller. However, with a sampling time of 5 ms, the system is much less stable: the loss of phase margin is far too great in this case.
Sampling effect:
2.1.3.2. Discrete-time PID Controller Implementation#
Since computers and microcontrollers cannot deal with continuous time, the control law has to be discretized. We’ll use
the parrallel form of a PID controller:
\(C(p)=\frac{O(p)}{E(p)}=K_p+\frac{K_i}{p}+K_dp\)
PID schema block:
Note: Other forms of PID controller are:
serial form \(C(p)=K_p(1+\frac{1}{\tau_ip})(1+\tau_sp)\)
standard form \(C(p)=K_p(1+\frac{1}{\tau_ip}+\tau_sp)\)
The time expression of this corrector is as follows:
\(o(t)=K_p.e(t) + K_i \int e(t)dt + K_d \frac{de(t)}{dt}\)
with the error: \(e(t)=SetPoint(t)-Input(t)\)
One of simplest methods available is to approximate this continuous-time ordinary differential equation by discrete-time difference equations or recurrence relations.
A first approach with parrallel form of the PID controller is:
for proportionnal term: \(K_p.e(t) \approx K_p.e_k\) where \(e_k = e(k.T_s)\)
for integration term: \(K_i\int e(t)dt \approx K_i\sum_{j=0}^{k} e_{j} T_s= K_i.I_k\) where \(I_k=I_k-1+e_{k}T_s\)
for derivation term: \(K_d\frac{de(t)}{dt} \approx K_d\frac{e_k - e_{k-1}}{T_s}\)
with the time step or sampling interval \(T_s\).
The PID output in response to measurements is thus:
\(I_k=I_{k-1}+e_{k}T_s\)
\(O_{k} = K_p e_{k} + K_i I_k + K_d\frac{e_k - e_{k-1}}{T_s}\)
Another more general approach applicable to any transfer function and therefore any form of corrector is as follows:
Write the transfer function in numerator/denominator form
\(C(p)=\frac{O(p)}{E(p)}=\frac{K_p p+K_i+K_d p^2}{p}\)Reformulate in the form of a differential equation
\(pO(p)=(K_p p+K_i+K_d p^2)E(p)\)
\(\dot o(t)=K_i e(t) + K_p \dot e(t) +K_d \ddot e(t)\)Discretize in the form of a recurrence equation
\(\frac{o_k-o_{k-1}}{T_q}=K_i e_k + K_p \frac{e_k-e_{k-1}}{T_s} + K_d \frac{e_k+ e_{k-2}-2e_{k-1}}{T_s^2}\)
\(o_k = o_{k-1} + K_i e_k T_s + K_p (e_k-e_{k-1}) + \frac{K_d}{T_s}(e_k+ e_{k-2}-2e_{k-1})\)
with
\(f_k=f(kT_s)\)
\(\dot f(t) \approx \frac{f(t)-f(t-T_s)}{T_s}\)
\(\ddot f(t) \approx \frac{f(t)+f(t-2T_s)-2f(t-T_s)}{T_s^2}\)
Note:
More advanced transformation or discretization method are avalaible.
More advanced PID controller forms can be find here
The following figure shows the effects of the digital implementation of the controller. Overall, the digital controller makes it possible to recover the behaviour of the continuous controller:
\(C(p)=\frac{O(p)}{E(p)}=K_p+\frac{K_i}{p}\)
using a simple recursive equation:
\(o_k = o_{k-1} + K_i e_k T_s + K_p (e_k-e_{k-1})\)
Digital controller:
2.1.4. Homework#
Install the Legacy Arduino IDE (1.8.X) + the Anaconda environment on your PC
Read this introduction to C code on Arduino.
2.1.5. References#
[Mihalič, 2022 ] Mihalič, F., Truntič, M., & Hren, A. (2022). Hardware-in-the-loop simulations: A historical overview of engineering challenges. Electronics, 11(15), 2462. Link
[Vanthuyne, 2009 ] Vanthuyne, T. (2009, September). An electrical thrust vector control system for the VEGA launcher. In Proceedings of the 13th European Space Mechanisms and Tribology Symposium (Vol. 670). Link
[Carpentier, 2011 ] Carpentier, B., Puel, F., Cros, C., Pin, B., Gerbe, L., & Solide, S. S. P. (2011). Stability analysis using a coupled model of TVC actuator control loop and structural dynamics of a Launch Vehicle. In 4th European Conference For Aerospace Sciences (EUCASS).
[Ellis, 2012] Ellis, G. (2012). Control system design guide: using your computer to understand and diagnose feedback controllers. Butterworth-Heinemann.