H-infinity Control of DC Motor Part I : Plant Modeling and Controller Synthesis

H_\infty 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 H_\infty 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, H_\infty 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 from Microchip as MCU.



Figure 1 a simple DC motor experimental setup

Step 1: 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

-->Pzsys_id
 Pzsys_id  = 
                                         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

-->P_id
 P_id  =
                                       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

-->roots(P_id.num)
 ans  =
    3436.017                
  - 240.67648 + 307.73657i  
  - 240.67648 - 307.73657i  
    200.                    
  - 94.567972

and poles at

-->roots(P_id.den)
 ans  =
  - 306.919 + 216.09937i  
  - 306.919 - 216.09937i  
  - 123.56295             
  - 8.4363                
  - 0.0078361

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

Note: 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.

Step 2: 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 |S(j\omega)| )

With some trial and error on weighting function adjustments, we are eventually satisfied with

(1)   \begin{equation*}  w_1(s) = \frac{0.5s+5}{s+0.01} \end{equation*}

(2)   \begin{equation*}  w_2(s) = \frac{0.01s+1}{0.004s+2} \end{equation*}

as weighting functions on S and T, respectively. In Scilab, these can be created as follows.

s=poly(0,'s');
a=0.002; 
m=2; 
wb=5;
w1_n=((1/m)*s+wb);
w1_d=(s+wb*a);
w1 = syslin('c',w1_n,w1_d);
 
w2_d=2*(s/500 + 1);
w2_n=s/100+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 S and T, The inverses of weighting functions (1) and (2) are plotted

gainplot([1/w1; 1/w2],0.001,100,['$1/W_1$';'$1/W_2$']);

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; 
        Bw1 zeros(size(Bw1,1),size(Bp,2));
        zeros(size(Bw2,1),size(Bw1,2)) zeros(size(Bw2,1),size(Bp,2))];
 
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
Pgen=syslin('c',Agp,Bgp,Cgp,Dgp);

Step 3: run the synthesis

A couple of Scilab commands could be used for H_\infty synthesis. hinf and h_inf, for example, are two separate commands with different syntax. To use the former,

[Ak,Bk,Ck,Dk]=hinf(Agp,Bgp,Cgp,Dgp,1,1,5);

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).

Step 4: reduce controller order

H_\infty 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

exec('kreduced.sci',-1);

Pack the controller if necessary

Sk=syslin('c',Ak,Bk,Ck,Dk);

and pass it to kreduced. You will see on the console some message such as below

-->Skr = kreduced(Sk);
Warning :
matrix is close to singular or badly scaled. rcond =    1.9157D-10
    150422.14  
    7713.887   
    115.87572  
    0.0339754  
    0.0003412  
    0.0000053  
    2.023D-09  
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 command

-->Skr_tf = ss2tf(Skr)
 Skr_tf  =
 
                                       2      3  
    384.79556 + 46179.923s + 1146.8162s - 500s   
    ------------------------------------------   
                                        2   3    
     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.

Step 5: 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 L, S, and T transfer functions from the plant model and the controller *after model reduction*

-->L=Skr_tf*Pqv_a; 
-->S=1/(1+L); 
-->T = 1-S;

and check the poles of T (or S)

-->roots(T.den)
 ans  =
 
  - 306.74046 + 215.68462i  
  - 306.74046 - 215.68462i  
  - 124.3588                
  - 13.073325 + 5.6429085i  
  - 13.073325 - 5.6429085i  
  - 6.9173679 + 1.1537611i  
  - 6.9173679 - 1.1537611i  
  - 0.0083338

The closed-loop frequency responses can be observed

-->gainplot([S; T],0.001,100,['$S(j\omega)$';'$T(j\omega)$']);

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 S and T

In H_\infty Control of DC Motor Part II : Controller Implementation we will discuss conversion of the controller to discrete-time and implementation on the target board.

Supplement

Scilab files used in this article

Comments

comments

Comments are closed.