is a recent linear control scheme that often sounds intimidating to a new control engineering student. And she might have in her mind already “Geez, this beast must require expensive control design software and run only on state-of-the-art laboratory hardware beyond my reach.” That is a total misconception indeed. I have been using Scilab functions related to robust control for some time and found them professional enough for real applications. I can analyze, synthesize, plot, simulate, and do whatever necessary with this open-source software to make sure that the resulting controller has the desired stability and performance properties. After satisfied, I implement my controller on a hand-wired prototype board that costs less than 50 bucks total. In this 2-part article, I am sharing my experience with you.
While I would not argue that the theory side of this so-called post-modern control could be mathematically involved, once one put some effort to grab the essence, the practical aspect is quite systematic and hence suitable for computer-aided analysis and design. In fact, the word “design” may be replaced by “synthesis” to indicate problem solvability with less human intervention. Of course, it is still mandatory for us to specify the requirements to a synthesis algorithm, typically in some form of weighting functions. The process afterwards is automatic. At the end we get our controller, or software complaint that none could meet our specs. Perhaps we have to relax this and that a little, and rerun the algorithm. The chance of success depends on plant complexity as well as how stringent the desired stability and performance criteria.
This article and the sequel address how to apply
control paradigm to a simple setup, angle control of a DC motor. The complete process is demonstrated, from constructing weighting functions from stability and performance requirements, forming generalized plant, synthesizing a controller, performing order reduction, and implementing on an embedded target board. Being in the model-based control category,
control requires a math model of the plant. So the first step is to acquire one using some type of system identification such as least-square method.
All the details involved were discussed in our previous articles. So this one is more like a how-to demonstration. For more insight or background, the following posts might be helpful
The DC motor is the same setup used in the DC Motor Speed Control Part I: Open-Loop Command
article, which is shown once again in Figure 1. The plant consists of those parts mounted on pine wood: DC motor, shaft encoder, and H-bridge driver. The controller is implemented on the target board on the right, using PIC24EP256MC202
Figure 1 a simple DC motor experimental setup
estimate a plant model
A function to generate PRBS input is implemented on the target board using (a C code example is provided on page 12 our System Modeling 101
article, then the input and output data is collected as shown in Figure 2. Note that the motor angle fluctuates around its initial angle at 180 degree by the PRBS input with strength +/- 5000. With 16-bit PWM and one direction bit, this translates to about 7% duty cycle, just enough to overcome shaft friction. Too large PRBS input strength might cause nonlinear motion or saturation.
Figure 2 PRBS input and motor angle response
Passing this input and output data to LSID algorithm in Scilab yields the following discrete-time plant
2 3 4
- 0.0000014 - 0.0000006z + 0.0000054z + 0.0000158z + 0.0000102z
2 3 4 5
- 0.0415585 + 0.1184330z + 0.2822806z + 0.1315280z - 1.4906741z + z
A standard controller synthesis is performed in continuous-time domain. Using a customized function as described in our Convert Transfer Function from Discrete to Continuous Time
article results in a continuous-time plant model
2 3 4 5
- 0.0009402 - 0.0000079s + 2.941D-08s + 1.146D-10s + 2.901D-13s - 9.479D-17s
2 3 4 5
- 0.0002877 - 0.0367550s - 0.0048099s - 0.0000557s - 0.0000002s - 2.500D-10s
with zeros at
- 240.67648 + 307.73657i
- 240.67648 - 307.73657i
and poles at
- 306.919 + 216.09937i
- 306.919 - 216.09937i
Note that the plant model is non-minimum phase, and has a low frequency pole in place of an integrator that it should have in theory. Its Bode plot is shown in Figure 3.
Figure 3 Bode plot of DC motor model obtained from LSID algorithm
All the process above is contained in Scilab script dcm_lsid.sce
Before we move on, it is worth mentioning that in a real application, one should verify that the model achieved is a good candidate of the physical plant. A way to do so is to feed another set of rich input data (different from the one used in the identification process) to both the plant and model, and then compare the output data to see how well they match. Since system identification is not the main focus of our discussion, I choose to omit this verification step to save space and time.
create weighting functions and form generalized plant
Standard S/T mixed sensitivity scheme is used for this problem. For initial design, let us make a crude stability and performance requirements, say
normalized tracking error less than 0.1 degree at low frequency region below 0.1 Hz
closed-loop bandwidth about 10 Hz
sufficient stability margin (no high peak on )
With some trial and error on weighting function adjustments, we are eventually satisfied with
as weighting functions on
, respectively. In Scilab, these can be created as follows.
w1 = syslin('c',w1_n,w1_d);
w2_d=2*(s/500 + 1);
w2 = syslin('c',w2_n,w2_d);
See H-infinity Synthesis with Scilab Part I : Problem Setup
article on the roles of parameters a, m, and wb and how to select them.
To observe the criteria on
, The inverses of weighting functions (1) and (2) are plotted
as shown in Figure 4. You can verify that they are fair candidates to the requirements listed above.
Figure 4 inverse of weighting functions
I cannot emphasize enough that weight selection is the crucial step of H_\infty control design. It is the last moment that we could tweak things to our desire, before the machine takes control. Some control theory background is helpful. If you are already good at classical control design, that same knowledge is inherited in the synthesis process by means of weighting functions. Otherwise, you may want to have a peek on Module 2
and Module 3
of our Scilab control engineering basics study modules.
Now that we have the plant model and weighting functions, a generalized plant can be formed. From my experience, converting all transfer functions to state-space systems and performing matrix augmentation gives best numerical result.
[Ap,Bp,Cp,Dp]=abcd(Pqv_a); // plant model (called P_id previously)
[Aw1,Bw1,Cw1,Dw1]=abcd(w1); // weighting function for S
[Aw2,Bw2,Cw2,Dw2]=abcd(w2); // weighting function for T
Agp=[Ap zeros(size(Ap,1),size(Aw1,2)) zeros(size(Ap,1),size(Aw2,2));
-Bw1*Cp Aw1 zeros(size(Aw1,1),size(Aw2,2));
Bw2*Cp zeros(size(Aw2,1),size(Aw1,2)) Aw2];
Bgp = [zeros(size(Bp,1),size(Bw1,2)) Bp;
Cgp = [-Dw1*Cp Cw1 zeros(size(Cw1,1),size(Cw2,2));
Dw2*Cp zeros(size(Cw2,1),size(Cw1,2)) Cw2;
-Cp zeros(size(Cp,1),size(Cw1,2)) zeros(size(Cp,1),size(Cw2,2))];
Dgp = [Dw1 0.001; 0 0;1 0]; // makes D12 full rank
run the synthesis
A couple of Scilab commands could be used for
, for example, are two separate commands with different syntax. To use the former,
where 1,1 are the number of control input and measured output. In this case our controller is SISO. The last argument is gamma variable you can adjust. Roughly speaking, smaller gamma means more stringent criteria. If you get a controller that yields unstable closed-loop, try relaxing gamma to larger value.
Note there is no iteration in hinf
. The command just runs once and returns a controller. In contrast,
[Sk,ro]=h_inf(Pgen, [1,1],0, 200000, 100);
iterates to yield a sub-optimal controller. Type help h_inf
to learn what the arguments are.
Regardless of which command used, at the end you get a controller in state-space format (packed or unpacked).
reduce controller order
synthesis generally returns a controller with higher order than necessary. From this DC motor example we get a controller of order 7. In an embedded application the resources are limited. High order means computation complexity and hence longer sampling period. Worse, the coefficients of higher order controller are subjected to numerical error, especially when using fixed-point, low performance MCUs.
So, we want to trim our controller to a lowest order one possible, without losing its properties significantly. How can we achieve such task?
A common approach is known as balanced truncation. Google it to learn the details. Here I just provide a Scilab function kreduced.sci
ready for use. Load this function to Scilab workspace
Pack the controller if necessary
and pass it to kreduced
. You will see on the console some message such as below
-->Skr = kreduced(Sk);
matrix is close to singular or badly scaled. rcond = 1.9157D-10
Enter number of states to keep:
The computed Hankel norms might be different depending on your controller. The idea is to find the first large gap in the values and truncate there. In the numbers above, we observe that between 115.87572 and 0.0339754, the controller states have noticable energy drop. This says the balanced states above 3 are quite insignificant. Hence we can reduce the controller order to 3. Type the number and press [Enter] to yield the reduced controller contained in Skr. To observe this controller in transfer function form, use Scilab ss2tf
-->Skr_tf = ss2tf(Skr)
384.79556 + 46179.923s + 1146.8162s - 500s
4.9087826 + 461.63448s + 31.256351s + s
Before getting into implementation on Part II, we end this article with the last step that should never be omitted.
checking closed-loop stability
An unstable feedback system is not only useless, it could also be hazardous to the environment and humans in the workplace. Despite what the theory says, synthesis algorithms may give a controller that is not stabilizing, mostly when the weighting functions are not chosen properly, or too stringent criteria requested. Or, a full-order controller from synthesis could be degraded by reduction scheme to the point that closed-loop stability is lost. In any case, we always have to check this property right after synthesis.
The easiest method for a continuous-time SISO feedback system is to check whether the closed-loop poles are Hurwitz; i.e., have negative real parts.
In Scilab, construct the
transfer functions from the plant model and the controller *after model reduction*
-->T = 1-S;
and check the poles of
- 306.74046 + 215.68462i
- 306.74046 - 215.68462i
- 13.073325 + 5.6429085i
- 13.073325 - 5.6429085i
- 6.9173679 + 1.1537611i
- 6.9173679 - 1.1537611i
The closed-loop frequency responses can be observed
as shown in Figure 5. Verify whether they meet our specs given earlier. If not, adjust the weighting functions and rerun the synthesis.
Figure 5 frequency responses of and
In Control of DC Motor Part II : Controller Implementation
we will discuss conversion of the controller to discrete-time and implementation on the target board.
Scilab files used in this article