The first step is to obtain the matrix A(q) characteristic polynomial A\left ( q\right ) =\begin{bmatrix} 0 & 1 & 0\\ 0 & 0 & 1\\ -\left ( 0.5+\varepsilon +3q_{1}+3q_{2}+2q_{1}q_{2}\right ) & -\left ( 1+q_{1}+q_{2}\right ) & -\left ( 1+q_{1}+q_{2}\right ) \end{bmatrix}
The method of polynomial over-bounding was tried first, but it was inconclusive. The attempt is included in the appendix. A graphical method was then tried based on the zero exclusion principle using set value of polytope of polynomials, but that also was inconclusive as the polygon seen crossing the zero as the frequency increased. The result of this attempt is described in the appendix.
Analysis using Hurwitz matrix Using P\left ( s\right ) =\left ( \varepsilon +3q_{1}+3q_{2}+2q_{1}q_{2}+\frac{1}{2}\right ) +s\left ( 1+q_{1}+q_{2}\right ) +s^{2}\left ( 1+q_{1}+q_{2}\right ) +s^{3} with 0\leq q_{1}\leq 1,0\leq q_{2}\leq 1. We setup the Hurwitz matrix for the above polynomial to find the conditions under which it is stable. H=\begin{bmatrix} 1+q_{1}+q_{2} & 1 & 0\\ 0.5+\varepsilon +3q_{1}+3q_{2}+2q_{1}q_{2} & 1+q_{1}+q_{2} & 0\\ 0 & 1+q_{1}+q_{2} & 1 \end{bmatrix}
Hence we only need to examine two cases. For \Delta _{1}=1+q_{1}+q_{2}, we see this is positive for all q, since 0\leq q_{i}\leq 1.
For \Delta _{2}=0.5-\varepsilon -q_{1}+q_{1}^{2}-q_{2}+q_{2}^{2} we need to determine the conditions which makes this positive. 0.5-\varepsilon -q_{1}+q_{1}^{2}-q_{2}+q_{2}^{2}>0
There are 20 increments, and for each \varepsilon _{i} an algebraic conditions was found on using the computer on q_{1},q_{2} which insures that \Delta _{2}>0. The result is tabulated below. In addition, 3D plot of \Delta _{2} is given, using q_{1},q_{2} as x,y and using the value of \Delta _{2} as the z-axis. This gives a visual view of the \Delta _{2} showing that it is indeed positive all the time using the constraints found on q_{i} for some specific \varepsilon used. A typical 3D for some \varepsilon is shown below for illustration
The above shows that for \varepsilon =0.005, \Delta _{2}>0, and hence the system is stable under the conditions 0\leq q_{1}\leq 0.4294,0\leq q_{2}\leq 1.
Robust stability depends on positive \Delta _{2}. In the above, we obtained conditions on \varepsilon ,q_{1},q_{2} which when met, will insure \Delta _{2}>0 and hence a stable polynomial. To visualize the solution in 3D, where we will use q_{1},q_{2},\varepsilon as the three axis, and use two colors: green to indicate the region where \Delta _{2}>0 and red color for the remaining region where \Delta _{2}\leq 0. Therefore, the space is a 3D cube enclosed in 1,1,\varepsilon . Using the above inequalities, the 3D plots are made.
Another 3D plot was made, using different 3D plot, called the surface plot, to help visualize the regions in different way. (some parts of the cube are not fully shown due to limited sampling used in producing the data).
We see from the above, that stability emerges for the region enclosed by 0\leq q_{1}\leq 0.5 and 0\leq q_{2}\leq 1 for low values of \varepsilon and for the region 0\leq q_{1}\leq 0.2 and 0\leq q_{2}\leq 1 for the large values of \varepsilon . The small \varepsilon is, and the smaller q_{1} is, the larger the stability region becomes.
Attempts were first made to obtain answer using zero exclusion condition. But both attempts resulted in inclusive result, as the polygon crosses the zero point. Algebraic reduction was used to obtain the constraints on solving for conditions on q_{1},q_{2} for making \Delta _{2}=0.5-\varepsilon -q_{1}+q_{1}^{2}-q_{2}+q_{2}^{2} positive as \varepsilon was incremented by small amount as shown above. Writing \Delta _{2} as \Delta _{2}=0.5-\varepsilon -\overbrace{\left ( q_{1}+q_{2}\right ) }^A+\overbrace{\left ( q_{1}^{2}+q_{2}^{2}\right ) }^B
When both q are around mid point of their range, for example q_{i}=\frac{1}{2}, then -A+B=-\frac{1}{2} and the result is \Delta _{2}=-\varepsilon , hence not stable.
As the q_{i} becomes larger than \frac{1}{2} then A does not dominate over B as much, but they are both larger now and the difference between A,B becomes smaller again as when they were both below \frac{1}{2}, which means now \Delta _{2} becomes larger, hence stable.
Here is a small table showing this variation. The above implies that the critical condition is where both q_{i} are close to each other in value. As they get closer to 0.5 then \varepsilon has to become smaller in order to keep \Delta _{2} positive.
Solving as constraint minimization problem Another attempt at theoretical support for the numerical computation, is to view this as constraint minimization problem, where we want to find the minimum of \Delta _{2} subject to constrain 0\leq q_{1}\leq 1,0\leq q_{2}\leq 1 and 0<\varepsilon <0.1, then the method of Lagrangian multipliers can be used.
Let the objective function be Q=\Delta _{2}=0.5-\varepsilon -q_{1}+q_{1}^{2}-q_{2}+q_{2}^{2} subject to 0\leq q_{1}\leq 1,0\leq q_{2}\leq 1,0<\varepsilon <0.1. Hence the Lagrangian L is\begin{align*} L & =Q+\lambda _{1}\left ( 1-q_{1}\right ) +\lambda _{2}\left ( 1-q_{2}\right ) +\lambda _{3}\left ( 0.1-\varepsilon \right ) +\lambda _{4}q_{1}+\lambda _{5}q_{2}+\lambda _{6}\varepsilon \\ & =0.5-\varepsilon -q_{1}+q_{1}^{2}-q_{2}+q_{2}^{2}+\lambda _{1}\left ( 1-q_{1}\right ) +\lambda _{2}\left ( 1-q_{2}\right ) +\lambda _{3}\left ( 0.1-\varepsilon \right ) +\lambda _{4}q_{1}+\lambda _{5}q_{2}+\lambda _{6}\varepsilon \end{align*}
Therefore \begin{align*} \frac{\partial L}{\partial q_{1}} & =-1+2q_{1}-\lambda _{1}+\lambda _{4}=0\\ \frac{\partial L}{\partial q_{2}} & =-1+2q_{2}-\lambda _{2}+\lambda _{5}=0\\ \frac{\partial L}{\partial \varepsilon } & =-1-\lambda _{3}+\lambda _{6}=0\\ \frac{\partial L}{\partial \lambda _{1}} & =1-q_{1}=0\\ \frac{\partial L}{\partial \lambda _{2}} & =1-q_{2}=0\\ \frac{\partial L}{\partial \lambda _{3}} & =0.1-\varepsilon =0\\ \frac{\partial L}{\partial \lambda _{4}} & =q_{1}=0\\ \frac{\partial L}{\partial \lambda _{5}} & =q_{2}=0\\ \frac{\partial L}{\partial \lambda _{6}} & =\varepsilon =0 \end{align*}
And \begin{bmatrix} 2 & 0 & 0 & -1 & 0 & 0 & 1 & 0 & 0\\ 0 & 2 & 0 & 0 & -1 & 0 & 0 & 1 & 0\\ 0 & 0 & 0 & 0 & 0 & -1 & 0 & 0 & 1\\ -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & -1 & 0 & 0 & 0 & 0 & 0 & 0\\ 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}\begin{bmatrix} q_{1}\\ q_{2}\\ \varepsilon \\ \lambda _{1}\\ \lambda _{2}\\ \lambda _{3}\\ \lambda _{4}\\ \lambda _{5}\\ \lambda _{6}\end{bmatrix} =\begin{bmatrix} 1\\ 1\\ 1\\ -1\\ -1\\ -0.1\\ 0\\ 0\\ 0 \end{bmatrix}
Plotting P_{\varepsilon } We now treat q_{1},q_{2} as random variables. Stability is still decided by \Delta _{2}=0.5-\varepsilon -q_{1}+q_{1}^{2}-q_{2}+q_{2}^{2}
P_{\varepsilon } was drawn for \varepsilon =0\cdots 1 to see how it shows up before zooming in.
The above shows that for 0\leq \varepsilon \leq 0.1 the probability of robust stability is high. We can zoom in to that region
We see that at \varepsilon =0.1 there is about 68\% chance that the system will be stable and for \varepsilon =0 the probability the system is stable is 100\%.
Finding the formula for the probability P_{\varepsilon } f_{Z}\left ( z\right ) is the probability density function of the random variable Z. To obtain f_{Z}\left ( z\right ) , we use the following two definitions. For random variable Z=X+Y where X is random variable drawn from f_{X} distribution and Y is random variable drawn from f_{Y} distribution, then the random variable Z will be drawn from distribution made from the convolution of f_{X} with f_{Y} f_{Z}\left ( z\right ) =\int _{-\infty }^{\infty }f_{X}\left ( z-x\right ) f_{Y}\left ( x\right ) dx
Let Z_{1}= -X-Y=-\left ( X+Y\right ), hence f_{Z_{1}}\left ( z\right ) =-\int _{-\infty }^{\infty }f_{U}\left ( z-x\right ) f_{U}\left ( x\right ) dx
The first method to try is the method of polynomial over-bounding If this method says the polynomial is robustly stable, then we are done. However, if this method says the polynomial is not robustly stable, then it can still be stable. Hence the polynomial over-bounding method is called inconclusive, and we need to try other methods.
But we start with this method since it is simple to use to check. Using the method of over-bounding, we first need to determine the bounds on \bar{q}_{i}. In other words, we convert the polynomial in q to an interval polynomial in \bar{q} \bar{q}_{0}^{-}=\min _{q\in Q}a_{0}\left ( q\right ) =\min _{0\leq q\leq 1}\left ( \varepsilon +3q_{1}+3q_{2}+2q_{1}q_{2}+\frac{1}{2}\right )
q_{1} | q_{2} | \varepsilon +3q_{1}+3q_{2}+2q_{1}q_{2}+\frac{1}{2} |
0 | 0 | \varepsilon +0.5 |
0 | 1 | \varepsilon +3.5 |
1 | 0 | \varepsilon +3.5 |
1 | 1 | \varepsilon +3+3+2+\frac{1}{2}=\varepsilon +8.5 |
Hence \bar{q}_{0}^{-}=\varepsilon +0.5
q_{1} | q_{2} | 1+q_{1}+q_{2} |
0 | 0 | 1 |
0 | 1 | 2 |
1 | 0 | 2 |
1 | 1 | 3 |
Hence \bar{q}_{1}^{-}=1
We now construct the four Kharitonov polynomials to check for stability using Hurwitz matrix method\begin{align*} K_{1}\left ( s\right ) & =\left ( \varepsilon +0.5\right ) +s+3s^{2}+s^{3}\\ K_{2}\left ( s\right ) & =\left ( \varepsilon +8.5\right ) +3s+s^{2}+s^{3}\\ K_{3}\left ( s\right ) & =\left ( \varepsilon +8.5\right ) +s+s^{2}+s^{3}\\ K_{4}\left ( s\right ) & =\left ( \varepsilon +0.5\right ) +3s+3s^{2}+s^{3} \end{align*}
Hence H_{1}=\begin{bmatrix} 1 & 1 & 0\\ \left ( \varepsilon +0.5\right ) & 3 & 0\\ 0 & 1 & 1 \end{bmatrix} and for stability we require that \Delta _{1}=1>0, \Delta _{2}=3-\left ( \varepsilon +0.5\right ) =2.5-\varepsilon >0 and \Delta _{3}=2.5-\varepsilon >0. Since \,0<\varepsilon <0.1, then we see that \Delta _{i},i=1\cdots 3 are all positive. Now we consider K_{2}.
H_{2}=\begin{bmatrix} 3 & 1 & 0\\ \left ( \varepsilon +8.5\right ) & 1 & 0\\ 0 & 3 & 1 \end{bmatrix} and for stability we require that \Delta _{1}=3>0, \Delta _{2}=-5.5-\varepsilon . We see that \Delta _{2}<0 since 0<\varepsilon <0.1. Hence
.
Therefore we need to try a different method.
Graphical method based on zero exclusion using set value of polytope of polynomials The P\left ( s\right ) =\left ( \varepsilon +3q_{1}+3q_{2}+2q_{1}q_{2}+\frac{1}{2}\right ) +s\left ( 1+q_{1}+q_{2}\right ) +s^{2}\left ( 1+q_{1}+q_{2}\right ) +s^{3} with 0\leq q_{1}\leq 1,0\leq q_{2}\leq 1 is multilinear in q. To use the value set for polytope of polynomials and apply graphical zero exclusion method in the hope to confirm robust stability, we have to convert the polynomial to an over-bounding affine linear polynomial in q by introducing new q_{3}=q_{1}q_{2}4
The polynomial P\left ( s\right ) =\left ( \varepsilon +3q_{1}+3q_{2}+2q_{3}+\frac{1}{2}\right ) +s\left ( 1+q_{1}+q_{2}\right ) +s^{2}\left ( 1+q_{1}+q_{2}\right ) +s^{3} with 0\leq q_{1}\leq 1,0\leq q_{2}\leq 1,0\leq q_{3}\leq 1 and 0<\varepsilon <\varepsilon ^{+}. Hence the uncertainty bounding set Q has 8 extremes q^{1}=\left ( 0,0,0\right ) ,q^{2}=\left ( 0,0,1\right ) ,q^{3}=\left ( 0,1,0\right ) ,q^{4}=\left ( 0,1,1\right ) ,q^{5}=\left ( 1,0,0\right ) ,q^{6}=\left ( 1,0,1\right ) ,q^{7}=\left ( 1,1,0\right ) ,q^{8}=\left ( 1,1,1\right ) . The eight associated polynomials generated are\begin{align*} P\left ( s,q^{1}\right ) & =\left ( \varepsilon +0.5\right ) +s+s^{2}+s^{3}\\ P\left ( s,q^{2}\right ) & =\varepsilon +2.5+s+s^{2}+s^{3}\\ P\left ( s,q^{3}\right ) & =\left ( \varepsilon +3.5\right ) +2s+2s^{2}+s^{3}\\ P\left ( s,q^{4}\right ) & =\left ( \varepsilon +5.5\right ) +2s+2s^{2}+s^{3}\\ P\left ( s,q^{5}\right ) & =\left ( \varepsilon +3.5\right ) +2s+2s^{2}+s^{3}\\ P\left ( s,q^{6}\right ) & =\left ( \varepsilon +5.5\right ) +2s+2s^{2}+s^{3}\\ P\left ( s,q^{7}\right ) & =\left ( \varepsilon +6.5\right ) +3s+3s^{2}+s^{3}\\ P\left ( s,q^{8}\right ) & =\left ( \varepsilon +8.5\right ) +3s+3s^{2}+s^{3} \end{align*}
From the above the nodes will be determined only by the unique polynomials. We will use only six out of the eight above, they are\begin{align*} P\left ( s,q^{1}\right ) & =\left ( \varepsilon +0.5\right ) +s+s^{2}+s^{3}\\ P\left ( s,q^{2}\right ) & =\varepsilon +2.5+s+s^{2}+s^{3}\\ P\left ( s,q^{3}\right ) & =\left ( \varepsilon +3.5\right ) +2s+2s^{2}+s^{3}\\ P\left ( s,q^{4}\right ) & =\left ( \varepsilon +5.5\right ) +2s+2s^{2}+s^{3}\\ P\left ( s,q^{7}\right ) & =\left ( \varepsilon +6.5\right ) +3s+3s^{2}+s^{3}\\ P\left ( s,q^{8}\right ) & =\left ( \varepsilon +8.5\right ) +3s+3s^{2}+s^{3} \end{align*}
We need to check that we have at least one nominal stable polynomial in the family of polynomials. Using P\left ( s,q^{1}\right ) =0.5+s+s^{2}+s^{3} as the stable member, we check it for stability first: H=\begin{bmatrix} 1 & 1 & 0\\ 0.5 & 1 & 0\\ 0 & 1 & 1 \end{bmatrix}
And\begin{align*} P\left ( j\omega ,q^{2}\right ) & =3.5+2j\omega +2\left ( j\omega \right ) ^{2}+\left ( j\omega \right ) ^{3}\\ & =3.5+2j\omega -2\omega ^{2}-j\omega ^{3}\\ & =3.5-2\omega ^{2}+j\left ( 2\omega -\omega ^{3}\right ) \end{align*}
Polynomial P\left ( s,q^{3}\right ) is the same as P\left ( s,q^{2}\right ) , so corner q^{2} and q^{3} map to same point in complex plane.\begin{align*} P\left ( j\omega ,q^{4}\right ) & =8.5+3j\omega +3\left ( j\omega \right ) ^{2}+\left ( j\omega \right ) ^{3}\\ & =8.5+3j\omega -3\omega ^{2}-j\omega ^{3}\\ & =8.5-3\omega ^{2}+j\left ( 3\omega -\omega ^{3}\right ) \end{align*}
And\begin{align*} P\left ( j\omega ,q^{5}\right ) & =\left ( \varepsilon ^{+}+0.5\right ) +j\omega +\left ( j\omega \right ) ^{2}+\left ( j\omega \right ) ^{3}\\ & =\left ( \varepsilon ^{+}+0.5\right ) +j\omega -\omega ^{2}-j\omega ^{3}\\ & =\left ( \varepsilon ^{+}+0.5\right ) -\omega ^{2}+j\left ( \omega -\omega ^{3}\right ) \end{align*}
And\begin{align*} P\left ( j\omega ,q^{6}\right ) & =\left ( \varepsilon ^{+}+3.5\right ) +2j\omega +2\left ( j\omega \right ) ^{2}+\left ( j\omega \right ) ^{3}\\ & =\left ( \varepsilon ^{+}+3.5\right ) +2j\omega -2\omega ^{2}-j\omega ^{3}\\ & =\left ( \varepsilon ^{+}+3.5\right ) -2\omega ^{2}+j\left ( 2\omega -\omega ^{3}\right ) \end{align*}
Polynomial P\left ( s,q^{7}\right ) is the same as P\left ( s,q^{6}\right ) , so corner q^{6} and q^{7} map to same point in complex plane.\begin{align*} P\left ( j\omega ,q^{8}\right ) & =\left ( \varepsilon ^{+}+8.5\right ) +3j\omega +3\left ( j\omega \right ) ^{2}+\left ( j\omega \right ) ^{3}\\ & =\left ( \varepsilon ^{+}+8.5\right ) +3j\omega -3\omega ^{2}-j\omega ^{3}\\ & =\left ( \varepsilon ^{+}+8.5\right ) -3\omega ^{2}+j\left ( 3\omega -\omega ^{3}\right ) \end{align*}
This diagram illustrates the mapping5
To find the cut off frequency \omega _{c}, using \omega _{c}=1+\frac{\max \left \{ q_{0}^{+},q_{1}^{+},\cdots ,q_{n-1}^{+}\right \} }{q_{n}^{-1}}, where here we use the interval polynomial \bar{p}\left ( s,\bar{q}\right ) =\left [ \varepsilon +0.5,\varepsilon +8.5\right ] +\left [ 1,3\right ] s+\left [ 1,3\right ] s^{2}+\left [ 1,1\right ] s^{3} found above. This results in \omega _{c}=1+\frac{\max \left \{ \varepsilon +8.5,3,3\right \} }{1}=1+\varepsilon +8.5=9.5+\varepsilon =9.51
The polygons generated did include the origin as the frequency \omega was increased and different \varepsilon tried. Below is an example running a program written to implement this method