A two-mass-one-spring system

A  two-mass-one-spring system: description of the problem

In this section, a simple example is used to illustrate how IQC\(\beta\) can be used to estimate the upper bound of the L2-induced gain of a pair of signals in an uncertain system. Consider the controlled two-mass-one-spring system in Figure 2.1. This example is adopted from the manual of MATLAB LMI Control Toolbox. It was a benchmark problem proposed in [5]. The goal of the original problem is to design an output-feedback law \(u = G_{c} (s)x_{2}\) that adequately rejects the disturbances \(w_{1}\), \(w_{2}\), and guarantees stability for values of the stiffness parameter \(k\) ranging in a certain given range. Here the control design problem will not be taken into account. We will use a proposed controller and illustrate how to verify stability of the interconnected control system by using IQC\(\beta\) .

Figure 2.1: Two-mass-one-spring system

For \(m_{1} = m_{2} = 1\), the dynamics of the system can be described by a fourth order state-space model


where the state vector \(x\) consists of the position and the velocity of masses \(m_{1}\) and \(m_{2}\). The matrices \(A\), \(B_{1}\) and \(B_{2}\) are given as follows

$$ A=\begin{bmatrix} 0 & 0 & 1 & 0\\0 & 0 & 0 & 1\\ -k & k & 0 & 0\\ k & -k & 0 & 0\end{bmatrix},\quad B_{1}=\begin{bmatrix}0\\0\\1\\0\end{bmatrix},\quad B_{2}=\begin{bmatrix} 0 \\0 \\ 0 \\ -1\end{bmatrix}$$

Matrix \(A\) depends affinely on \(k\). By letting \(k = 1.25+\frac{3}{4}\bar{k}\), one can be further expressed as \(A = A_{0} +A_{1}\bar{k}A_{2}^{T}\), where \(A_{0}\), \(A_{1}\), \(A_{2}\) are again constant matrices

$$A_{0}=\begin{bmatrix} 0 & 0 & 1 & 0\\0 & 0 & 0 & 1\\ -1.25 & 1.25 & 0 & 0\\ 1.25 & -1.25 & 0 & 0\end{bmatrix},\quad A_{1}=\begin{bmatrix} 0\\0\\-\frac{\sqrt{2}}{2}\\\frac{\sqrt{2}}{2}\end{bmatrix},\quad A_{2}=\begin{bmatrix} \frac{3\sqrt{2}}{4}\\ -\frac{3\sqrt{2}}{4}\\0\\0\end{bmatrix},$$

and \(\bar{k}\) is an unknown constant in the range of \([−a,\ a]\). When \(a = 1\), the spring coefficient \(k\) belongs to the range of \([0.5,\ 2]\), which is the control design specification of the original problem in [5].

A fourth-order stabilizing controller \(u = C_{c}(sI-A_{c})^{-1}B_{c}x_{2}\) was proposed in [5], where

$$ A_{c}=\begin{bmatrix} 0 & -0.7195 & 1 & 0\\0 & -2.9732 & 0 & 1\\ -2.5133 & 4.8548 & -1.7287 & -0.9616\\ 1.0063 & -5.4097 & -0.0081 & 0.0304\end{bmatrix},\quad B_{c}=\begin{bmatrix} 0.720\\ 2.973\\ -3.37\\ 4.419\end{bmatrix},\quad C_{c}^{T}=\begin{bmatrix} -1.506\\ 0.494\\ -1.738\\ -0.932\end{bmatrix}$$

The mass-spring control system can be expressed as the block diagram shown in Figure 2.2. where \(B_{1}\), \(B_{2}\), \(C\) are \([0 \ 0 \ 1 \ 0]^{T}\), \([0 \ 0 \ 0 \ -1]^{T}\), \([0 \ 1 \ 0 \ 0]\), respectively.

 Figure 2.2: Block diagram of the two-mass-one-spring system

In the following, we will show how to estimate the gain from \([w1,\ w2]\) to \(x_{2}\) by IQC\(\beta\). This is done by computing an estimate of the energy gain from \([w1,\ w2]\) to \(x_{2}\). At current stage, checking stability is always performed by estimating the gain from some external disturbance to any internal signal in the system. A finite gain implies stability. We assume that matrices \(A_{0}\), \(A_{1}\), \(A_{2}\), \(B_{1}\), \(B_{2}\), and \(C\), as well as the LTI systems \(G_{p}(s)\) and \(G_{c}(s)\), are already defined in the MATLAB workspace. We will also assume that the readers have experience of working with MATLAB Control Systems Toolbox.

First Step - Initialize abstract IQC-environment

The first step is to initialize the abstract IQC-environment. This is done by typing

>> abst_init_iqc;

The M-file abst_init_iqc.m tells MATLAB the "rules" of the abstract IQC-environment: which types of abstract data can be defined, which type of external data is allowed, and which operations can be done in the abstract environment. These rules are to be discussed in details in the later chapters. Note that initialization needs to be done before using any command provided by IQC\(\beta\) .

Second Step - Describe the system

The next step is to describe the system, i.e. to define how the signals in the system are related. First of all, we need to identify the blocks which correspond to the \(\Delta_{k}\) blocks mentioned in the previous section. These blocks should be excluded at this stage, and later be replaced by set-valued functions defined using IQCs. In this example,\(\bar{k}\) is uncertain and hence the \(\bar{k}\) block is the one which should be removed. Secondly, we have to identify the so-called "basic signals" of the system. A basic signal is a signal that cannot be derived from any other signals through a constant LTI transformation. Usually, they are external disturbances and the outputs of uncertain blocks. In this example, \(w_{1}\), \(w_{2}\) and \(w_{3}\) are basic signals.

Sometimes, in order to facilitate the process of describing the system, extra basic signals are defined as dummies for signals which can be represented as LTI transformations of the other signals. Then, these signals are declared to be the same as their substitutes by using "==". For example, the two-mass-one-spring system in Figure 2.2 has a loop in which every signal is a LTI transformation of the others. One way to describe this loop is to break the loop as shown in Figure 2.3 and introduce an extra basic signal u. We can then define \(x\), \(x_{2}\), and \(u_{c}\) as LTI transformations of basic signals including \(u\). In the end, we declare that \(u\) is the same as \(u_{c}\) by using "==". Therefore, we need four basic signals in this example: \(w_{1}\), \(w_{2}\), \(w_{3}\), and \(u\).

Figure 2.3: Constant LTI part of the two-mass-one-spring system

To define these basic signals, one uses the M-function signal.m

>> w1=signal;

>> w2=signal;

>> w3=signal;

>> u=signal;

Here signal.m is the M-function used to declare basic signals. The syntax is

>> identifier_of_the_signal = signal(size_of_the_signal)

For example, one can define a vector signal w with size 4 by typing

>> w=signal(4);

The default size of the signalis 1.

After defining basic signals, we can go forward and relate the remaining signals (\(x\), \(x_{2}\), \(v\), and \(u_{c}\)) in the system to the basic signals we just defined.

>> x=Gp*(A1*w3+B1*(u+w1)+B2*w2);

>> v=A2*x;

>> x2=C*x;

>> uc=Gc*x2;

Note that \(x_{2}\) is the second element of \(x\), and therefore, we can also use subscripted reference

>> x2=x(2);

instead of making x2=C*x.

To close the control loop, we need to tell MATLAB that \(u\) and \(u_{c}\) are identical. This is done by the operation "=="

>> u==u_c;

The meaning of "==" is very different from the meaning of "=". If "=" would have been used instead, the result would simply be re-assigning \(u\) as an identifier to the signal identified by \(u_{c}\), while the signal originally identified by \(u\) still remains in the memory (and has no identifier now!). With "==", the command \(u==u_{c}\) tells MATLAB that \(u\) and \(u_{c}\) are identical. This operation is called linking and should be used with caution. Wrong usage of this function will result in error. For example, a command

>> u==2*u;

will most likely result in an error at the stage of running the optimizer iqc_gain_tbx.m.

Third Step - IQC

The next step is to relate \(w_{3}\) and \(v\) with a set-valued function which is defined using integral quadratic constraints. Since \(|\bar{k}| \leq a\), we know that

$$ \int_{0}^{\infty} (a^{2}|v(t)|^{2}-|w_{3}(t)|^{2})dt \geq 0.$$

This inequality remains valid after multiplication by any scalar larger than zero. Hence, the following parameterized integral inequality holds for all pairs of (\(w_{3}\), \(v\)) which satisfy \(w_{3} =\bar{k}v\)

$$ \int_{0}^{\infty} M \cdot (a^{2}\cdot v^{2} - w_{3}^{2})dt \geq 0.$$

where \(M \geq 0\) is the corresponding λ parameter mentioned in the previous section. To define this relationship, type

>> M=variable;

>> M>0;

>> v'*(aˆ2*M)*v-w3'*M*w3>0;

Here variable.m, which gives an analog of type 3 declaration in LMI Control Toolbox, is used to define the variable M.

The syntax of this M-function is

>> identifier_of_the_variable = variable(structure_of_the_variable)

where structure of the variable follows exactly the convention of the LMI Control Toolbox, and has the default value 1. Since M is a scalar variable, we simply use the default value. There are other functions which allow users to define special types matrix variables (for example, symmetric and skew symmetric matrices).

The description M>0 tells MATLAB that M can only be greater than (or equal to) 0, while v'*(aˆ2*M)*v-w3'*M*w3>0 describes exactly the IQC. Note that, for the sake of simplicity, "\(>\)"  and "\(<\)" are used to replace "\(\geq\)" and "\(\leq\)" in IQC\(\beta\).

Fourth Step -Estimate the L2-gain

The last step of the analysis is to execute the optimizer to estimate the L2-gain from [w1; w2] to x2. This is done by the M-function iqc_gain_tbx.m.

>> g=iqc_gain_tbx([w1;w2],x2)

The execution of iqc_gain_tbx trigger the following events: first, the IQC Toolbox scripts input between the commands abst_init_iqc and g=iqc_gain_tbx([w1;w2],x2) are collected and processed to form an optimization problem which corresponds to the worst-case L2-gain estimation problem. The optimization problem is in the form of the semi-definite program (SDP). Secondly, the genetic SDP solver coming with the MATLAB LMI Control Toolbox is called to solve the optimization problem. 

More precisely, an attempt will be made to choose the variable parameters (M and g in this case, where g is defined inside the M-function iqc_gain_tbx.m) so that

$$ \int_{0}^{\infty} (g(w_{1}^{2}+w_{2}^{2})-x_{2}^{2})dt > \int_{0}^{\infty} M(a^{2}v^{2}-w_{3}^{2})dt$$

holds for all non-zero L2 signals \(w_{1}\), \(w_{2}\), \(w_{3}\). In fact, the attempt to minimize g will also be made, and the SDP solver will try to find the smallest possible value of g.

If the system of LMIs will be found infeasible, the output g will be empty. Otherwise, an upper bound of g will be produced. The optimization will also produce some “optimal” value of M. This value can be displayed by typing

>> value_iqc(M)

Note that if the LMIs are found infeasible, execution of the command value_iqc(M) will result in an error message.

Full MATLAB Descriptions

Following is the summary of all MATLAB descriptions for this robustness problem: (we assume that \(A_{0}\), \(A_{1}\), \(A_{2}\), \(B_{1}\), \(B_{2}\), \(C\), \(G_{p}(s)\), and \(G_{c}(s)\) are already defined.)

>> abst_init_iqc;

% (1) initialize abst iqc environment

>> a=0.1;

% (2) set the value of a

>> w1=signal;

% (3) declare basic signal w1

>> w2=signal;

% (4) declare basic signal w2

>> w3=signal;

% (5) declare basic signal w3

>> u=signal;

% (6) declare basic signal u

>> x=Gp*(A1*w3+B1*(u+w1)+B2*w2);

% (7) define signal x

>> v=A2*x;

% (8) define signal v

>> x2=C*x;

% (9) define signal x2

>> uc=Gc*x2;

% (10) define signal uc

>> M=variable;

% (11) declare variable M

>> M>0;

% (12) impose a constraint on M

>> v’*(aˆ2*M)*v-w3’*M*w3>0;

% (13) define an IQC

>> u==uc;

% (14) declare u and uc are identical

>> g=iqc_gain_tbx([w1;w2],x2)

% (15) estimate the L2 gain

Here, the estimate of the L2-gain from [w1; w2] to x2 is 17.2881.

It is finite, and we conclude that the controlled system is stable for \(k\) ranging in \([0.5,\ 2]\).

Finally, to obtain the “optimal” value of multiplier M, we can use the following commands

>> iqc_value

>> value_iqc(M)

In this case, the result is 117.9644, which means the value of M to obtain L2-gain estimation 17.4074 is 117.9644.