Part I: Discretization and SimulationOne approach quite common in industrial control nowadays, when doing analysis and design with a real plant is expensive, is to implement the plant model on a computer, usually an embedded system. Some advantages of this so-called Hardware-In-the-Loop (HIL) approach, in addition to lowering development cost and time, are flexibility in adjusting the parameters, and the ability to perform extreme experiments which might damage the real plant.
From our previous articles on robot joint driven by DC motor with harmonic drive mechanism (from now I will abbreviate it as HDM for short) , we obtain a linear model and simulate with PID and custom-designed controllers. Here we continue our discussion on how to implement this math model on a real-time embedded system, to be used as an HIL plant for feedback control analysis and design. To be honest, this is my first HIL project that is based on a math model of a physical plant, except for some simple memoryless systems. Certain factors, both hardware and software related, such as ADC and PWM accuracy, delay, fixed-point math and other numerical issues, could affect the performance of this embedded system. With limited resources, at the end I am not so sure it could resemble the real plant good enough to warrant its place in a control lab. Nevertheless, it is worth a try. If things work out properly, it could be an interesting topic for my control engineering students for the upcoming semester. If not, I believe we learn something together in such attempt. 😉 The proposed procedure is as follows
- Discretize the continuous-time transfer function of the HDM
- Simulate the discrete-time plant and compare the responses with continuous-time model
- Adjust overall plant gain to suit the input and output limits of embedded system
- Implement on the selected microcontroller (in this case, a PIC24EP product from Microchip)
- Experiment with the HIL model
DiscretizationFrom the HDM block diagram derived from our previous article and displayed again in Figure 1, we can manipulate to achieve a reduced block diagram in Figure 2, consisting of only 2 sub-blocks
(5)with sampling period , we eventually have the discrete-time representation of and as
(8)Once again, the routine work is left as an exercise.
Verification with ScilabBefore doing anything else, it makes sense to verify that our derivation of the discrete-time transfer functions (6), (7) and their coefficients (8) are correct. Here we show how to verify (6) using discretization command in Scilab. You can use similar process for (7). Note :all Scilab commands listed below are included in the script file hdminit.sce Given a set of parameters, we construct the continuous-time plant as before.
Executing the above script yields the continuous-time transfer function
km = 100; // torque constant kb = 1; // back EMF constant k = 1000; // torsional stiffness of harmonic drive r = 10; // gear ratio L = 0.1 // armature inductance R = 1; // armature resistance Jm = 1; // motor inertia Bm = 0.01; // motor shaft friction Jl = 3; // load inertia Bl = 0.05; // load friction // plant denominator coefficient ac5 = L*Jm*Jl; ac4 = R*Jm*Jl+L*(Jm*Bl + Jl*Bm); ac3 = L*(Jm*k + Bm*Bl + Jl*k) + R*(Jm*Bl + Jl*Bm) + r*km*kb*Jl; ac2 = L*k*(Bm+Bl)+R*(Jm*k + Bm*Bl + Jl*k) + r*km*kb*Bl; ac1 = k*(L*k + R*(Bm+Bl) - k*L + r*km*kb); s=poly(0,'s'); pl = Jl*s^2+Bl*s+k; pm = Jm*s^2+Bm*s+k; // transfer function from V to load angle pnum1 = km*k; pden1 = ac5*s^5 + ac4*s^4 + ac3*s^3 + ac2*s^2 + ac1*s;; P1 = syslin('c',pnum1,pden1);
(9)Specify the sampling period of discrete-time system, say, 10 milliseconds.
To use the Scilab command cls2dls to discretize the plant, first it must be described in state-space form
T = 0.01;
Then put this plant and sampling time as arguments to cls2dls, which outputs a discrete-time state space system using Bilinear transformation.
P1ss = tf2ss(P1);
For comparison with our hand computation, we prefer a discrete-time transfer function. There is some other way to do this, but here we simply use the control textbook relationship .This is done in Scilab as follows
Pd1ss = cls2dls(P1ss, T);
Running the script in Scilab yields
[Ad,Bd,Cd,Dd]=abcd(Pd1ss); z = poly(0,'z'); Pd1 = Cd*inv(z*eye(size(Ad,1),size(Ad,2)) - Ad)*Bd + Dd;
Now we construct another discrete-time plant from our hand calculation.
-->Pd1 Pd1 = 2 3 4 5 0.0000008 + 0.0000039z + 0.0000078z + 0.0000078z + 0.0000039z + 0.0000008z -------------------------------------------------------------------------- 2 3 4 5 - 0.9224925 + 3.8350285z - 6.9849487z + 7.1298534z - 4.0574407z + z
To make a comparison, notice that the coefficient of highest in the denominator of Pd1 is normalized to 1, but in p1denz it equals alz5. So we have to normalize both p1numz and p1denz with this quantity. Using syslin command for this discrete-time plant yields the same transfer function as the one computed by Scilab cls2dls command.
// numerator coefficients of p1(z) b1z5 = T^5; b1z4 = 5*T^5; b1z3 = 10*T^5; b1z2 = 10*T^5; b1z1 = 5*T^5; b1z0 = T^5; // denominator coefficients of p1(z) a1z5 = 32*ac5 + 16*T*ac4 + 8*T^2*ac3 + 4*T^3*ac2 + 2*T^4*ac1; a1z4 = -160*ac5 - 48*T*ac4 - 8*T^2*ac3 + 4*T^3*ac2 + 6*T^4*ac1; a1z3 = 320*ac5 + 32*T*ac4 - 16*T^2*ac3 - 8*T^3*ac2 + 4*T^4*ac1; a1z2 = -320*ac5 + 32*T*ac4 + 16*T^2*ac3 - 8*T^3*ac2 - 4*T^4*ac1; a1z1 = 160*ac5 - 48*T*ac4 + 8*T^2*ac3 + 4*T^3*ac2 - 6*T^4*ac1; a1z0 = -32*ac5 + 16*T*ac4 - 8*T^2*ac3 + 4*T^3*ac2 - 2*T^4*ac1; // ----- p1(z) ------------------- p1numz = km*k*(b1z5*z^5 + b1z4*z^4 + b1z3*z^3 + b1z2*z^2 + b1z1*z + b1z0); p1denz = a1z5*z^5 + a1z4*z^4 + a1z3*z^3 + a1z2*z^2 +a1z1*z + a1z0;
and this should make us happy for the moment.
-->Pd1h = syslin('d',p1numz/a1z5,p1denz/a1z5) Pd1h = 2 3 4 5 0.0000008 + 0.0000039z + 0.0000078z + 0.0000078z + 0.0000039z + 0.0000008z -------------------------------------------------------------------------- 2 3 4 5 - 0.9224925 + 3.8350285z - 6.9849487z + 7.1298534z - 4.0574407z + z
Step Response SimulationTo see how well the resulting discrete-time transfer functions (6) and (7) represent their continuous-time counterparts, we must simulate the two and compare their responses. Before running any simulation model below, first you need to execute the script hdminit.sce to initialize all the variables used in the diagram. Figure 3 shows an Xcos diagram pid_lf_cdcompare.zcos, a PID control comparison with load angle as feedback. The continuous time and discrete time plants are placed in the upper and lower feedback loops, respectively. So the upper loop is purely continuous system, while the lower loop is a hybrid system (continuous-time PID controller, discrete-time plant).
Adjusting Plant Gain to Embedded HardwareSo far we have developed a discrete-time model quite ready to be implemented on a microcontroller. There is one issue left to be considered, and is fundamentally related to the properties of hardware used. We ask a simple question. How can you be sure that the plant gain is comparable with the range of input and output of your hardware? The extreme cases are too high gain that easily causes output saturation, or too small gain that requires too excessive input to drive the system than the controller could deliver. So we need to calibrate the plant gain to take advantage of the full input and output ranges. For the plant transfer functions described by (6),(7), the plant gain is a complicated function of many parameters. In fact, due to the presence of integrator the DC gain is infinity. So if we simulate an open-loop step response, the motor does not settle at a value but continues to rotate. This makes gain adjustment somewhat tricky. Nevertheless, it is not too difficult when we realize that, in a real motor, a property that is already available as standard specs is the velocity, commonly specified in Revolution Per Minute (RPM). In a simulation diagram, motor velocity can be computed by adding a derivative block as shown in Figure 7. Note in the diagram that we “measure” both and and their derivatives. It is however only one variable, in this case , the motor shaft velocity, we want to calibrate, since two angle outputs are proportionally related by the gear ratio.
(10)Note that the term in parenthesis represents the transfer function from input voltage to motor shaft velocity. Let represent output rate of change per second, then at DC we have
(11)At maximum input, we want the motor to rotate at maximum speed, and translate this to maximum output. We also have to take into consideration the motor sense, clockwise or counter-clockwise. So the range for each direction is half of the input range, say, for ADC input value of 511, the motor rotates at RPM, and the PWM maximum output is 4095. Substituting all these quantities and do some math, we have
(12)For example, a motor with maximum speed 3,000 RPM yields . With this gain and step input of 511, the output changes at rate 204,350 units/sec. The equation (12) is coded in hdminit.sce and used in Xcos model speedtuning.zcos. Figure 9 shows the simulation result of the output rate, when the step input is applied at t = 1 second. We see the fluctuation around the expected value.
Scilab Script and Model FilesThe files used this document are
- hdminit.sce : Scilab script file for parameter initialization and transfer function computation
- pid_mf_cdcompare.zcos : Xcos model comparison for PID with motor angle feedback
- pid_lf_cdcompare.zcos : Xcos model comparison for PID with load angle feedback
- speedtuning.zcos : Xcos model for plant gain adjustment
Note :execute the script file first before running simulationDownload all as zip file: hilhdm1.zip
- Show derivation details of continuous-time transfer functions and their coefficients (1) – (4).
- Show derivation details of discrete-time transfer functions and their coefficients (6) – (8).
- Derive (12).
- Experiment with Xcos model comparisons. Try changing sampling time, and PID gains. Can you explain their effects on the responses?