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} \] Therefore \begin{align*} p\left ( s,q\right ) & =\left \vert sI-A\left ( q\right ) \right \vert \\ & =\begin{vmatrix} s & -1 & 0\\ 0 & s & -1\\ \left ( 0.5+\varepsilon +3q_{1}+3q_{2}+2q_{1}q_{2}\right ) & \left ( 1+q_{1}+q_{2}\right ) & s+\left ( 1+q_{1}+q_{2}\right ) \end{vmatrix} \\ & =\left ( \varepsilon +3q_{1}+3q_{2}+2q_{1}q_{2}+\frac{1}{2}\right ) +\left ( 1+q_{1}+q_{2}\right ) s+\left ( 1+q_{1}+q_{2}\right ) s^{2}+s^{3}\\ & =a_{0}\left ( q\right ) +a_{1}\left ( q\right ) s+a_{2}\left ( q\right ) s^{2}+a_{3}\left ( q\right ) s^{3} \end{align*}
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} \] The leading minors are\begin{align*} \Delta _{1} & =1+q_{1}+q_{2}\\ \Delta _{2} & =0.5-\varepsilon -q_{1}+q_{1}^{2}-q_{2}+q_{2}^{2}\\ \Delta _{3} & =\Delta _{2} \end{align*}
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 \] In other words, the minimum of \(\Delta _{2}\) should be positive to insure stability. This is global minimization with constrain problem. However, an algebraic reduce method was used instead to obtain the limits on \(q_{1},q_{2}\) for each different \(\varepsilon _{i}\) where \(\varepsilon _{i}\) was incremented by \(0.005\) from \(0\) to \(0.1\)
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 \] Since each \(0\leq q_{i}\leq 1\) then \(q_{i}^{2}<q_{i}\). This means that the smaller the \(q_{i}\) becomes then \(A\) will dominate over \(B\) more in size but they are both small, and hence \(\Delta _{2}\) is positive .
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} \] Solving the above \(Ax=b\) system using least squares gave the solution as \(q_{1}=0.5,q_{2}=0.5,\varepsilon =0.05\), which is where the minimum is. Using the solution we can find \(\Delta _{2}\) at these values \[ \Delta _{2}=-0.05 \] The above Lagrangian multiplier method to find the minimum was implemented in separate program allowing one to change the limits of \(\varepsilon \) and \(q_{1},q_{2}\) to see where the minimum shows up for each different combination. Here are few screen shots showing different results. This method was used to verify the numerical 3D based plots shown above by verifying the stable points are where they are shown in the 3D plots.
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}\] If we call \(\Delta _{2}\) as the random variable \(Z\) which is now function of random variables \(q_{1},q_{2}\) renamed to be \(X,Y\) and are drawn from uniform distribution, then we can write \[ Z=0.5-\varepsilon -X+X^{2}-Y+Y^{2}\] Then \[ \Pr \left ( z_{1}\leq Z\leq z_{2}\right ) ={\displaystyle \int \limits _{z_{1}}^{z_{1}}} f_{Z}\left ( z\right ) dz \] or\[ \Pr \left ( Z\leq z\right ) =F_{Z}\left ( z\right ) ={\displaystyle \int \limits _{-\infty }^{z}} f_{Z}\left ( z\right ) dz \] Where \(F_{Z}\left ( z\right )\) is the cumulative distribution function. To find \(\Pr \left ( Z>0\right )\) then\[ \Pr \left ( Z>0\right ) ={\displaystyle \int \limits _{0}^{\infty }} f_{Z}\left ( z\right ) dz \] \(\varepsilon \) is fixed each time before finding \(f_{z}\left ( z\right ) \). Hence for each \(\varepsilon \) there will be a different \(f_{Z}\left ( z\right )\) which we then use to find \(\Pr \left ( Z>0\right ) \) for robust stability, since \(\Pr \left ( Z>0\right ) \) is the same as asking what is the probability that \(\Delta _{2}>0\) for a given \(\varepsilon \).
\(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 \] Where \(f_{X}\) is the pdf of the uniform distribution for \(Q\) which is defined for \(q=\left [ 0,1\right ] \) as\[ f_{X}\left ( x\right ) =\left \{ \begin{array} [c]{cc}1 & 0\leq x\leq 1\\ 0 & \text{otherwise}\end{array} \right . \] For \(Z=XY\) evaluation, we need a product of two random variables. This is given by\[ f_{Z}\left ( z\right ) =\int _{-\infty }^{\infty }f\left ( x\right ) f\left ( z/x\right ) \frac{1}{\left \vert x\right \vert }dx \] Using the above, we can find that \(f_{Z}\left ( z\right ) \) for \(Z=0.5-\varepsilon -X+X^{2}-Y+Y^{2}\).
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 \] Now let \(Z_{2}=X^{2}\). Using the product formula, we write \[ f_{Z_{2}}\left ( z\right ) =\int _{-\infty }^{\infty }f_{U}\left ( x\right ) f_{U}\left ( z/x\right ) \frac{1}{\left \vert x\right \vert }dx \] Therefore for \(Z_{3}=-X-Y+X^{2}\) we now have \(Z_{3}=Z_{1}+Z_{2}\) and now we use the addition formula\[ f_{Z_{3}}\left ( z\right ) =\int _{-\infty }^{\infty }f_{Z_{1}}\left ( z-x\right ) f_{Z_{2}}\left ( x\right ) dx \] For \(-Y^{2}\), let \(Z_{4}=-Y^{2}\) and using the product formula gives\[ f_{Z_{4}}\left ( z\right ) =-\int _{-\infty }^{\infty }f_{U}\left ( x\right ) f_{U}\left ( z/x\right ) \frac{1}{\left \vert x\right \vert }dx \] Hence we now have the following \(Z_{5}=Z_{3}+Z_{4}\) and the pdf is \[ f_{Z_{5}}\left ( z\right ) =\int _{-\infty }^{\infty }f_{Z_{3}}\left ( z-x\right ) f_{Z_{4}}\left ( x\right ) dx \] Finally, we have \(Z=0.5-\varepsilon +Z_{5}\) which have the pdf\[ f_{Z}\left ( z\right ) =\frac{1}{\left \vert 0.5-\varepsilon \right \vert }f_{Z_{5}}\left ( z-\left ( 0.5-\varepsilon \right ) \right ) \] With the help of the computer, formula for the pdf of \(f_{Z}\left ( z\right ) \) was obtained which was used to generate the above plots of \(P_{\varepsilon }\)\[ f_{Z}\left ( z\right ) =\left \{ \begin{array} [c]{ccc}\pi & & 0<q+\epsilon \leq 0.25\\ 2\left ( \operatorname{arcsec}\left ( 2\sqrt{q+\varepsilon }\right ) -\arctan \left ( \sqrt{-1+4q+4\varepsilon }\right ) \right ) & & 0.25\leq q+\varepsilon <0.5\\ 0 & & \text{otherwise}\end{array} \right . \]
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 ) \] Setting up the following table
| \(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 \] And\[ \bar{q}_{0}^{+}=\max _{q\in Q}a_{0}\left ( q\right ) =\varepsilon +8.5 \] And\[ \bar{q}_{1}^{-}=\min _{q\in Q}a_{1}\left ( q\right ) =\min _{0\leq q\leq 1}\left ( 1+q_{1}+q_{2}\right ) \] Setting up a table
| \(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 \] And\[ \bar{q}_{1}^{+}=\max _{q\in Q}a_{1}\left ( q\right ) =3 \] And similarly, \(\bar{q}_{2}^{-}=1,\bar{q}_{2}^{+}=3\). And for \(\bar{q}_{3}\), since it has no uncertainties, then \(\bar{q}_{3}^{-}=\bar{q}_{3}^{+}=1\), Hence the over-bounding interval polynomial is\begin{align*} \bar{p}\left ( s,\bar{q}\right ) & =\left [ \bar{q}_{0}^{-},\bar{q}_{0}^{+}\right ] +\left [ \bar{q}_{1}^{-},\bar{q}_{1}^{+}\right ] s+\left [ \bar{q}_{2}^{-},\bar{q}_{2}^{+}\right ] s^{2}+\left [ \bar{q}_{3}^{-},\bar{q}_{3}^{+}\right ] s^{3}\\ & =\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} \end{align*}
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} \] For \(\Delta _{1}=1>0\), and \(\Delta _{2}=0.5>0\) and \(\Delta _{3}=\Delta _{2}>0\) as well. Hence we verified the stable member exist. Now we need to generate the polygonal value set for each of the polynomials above. \begin{align*} P\left ( j\omega ,q^{1}\right ) & =0.5+j\omega +\left ( j\omega \right ) ^{2}+\left ( j\omega \right ) ^{3}\\ & =0.5+j\omega -\omega ^{2}-j\omega ^{3}\\ & =0.5-\omega ^{2}+j\left ( \omega -\omega ^{3}\right ) \end{align*}
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 \] Therefore the sweep frequency from be \(0<\omega <9.51\) for the simulation6 . A program was written to simulate the value sets for this problem. For each \(\varepsilon \) value, the value set is displayed over the sweep frequency to see if the polygon will include the origin or not. \(\varepsilon \) was changed by small increments and the simulation was run again.
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
