Stability and robustness are fundamental design requirements of any control system. Consequently, stability analysis is a vital stage in the design and development process of a control system. In addition to providing information about the inherent stability of a system, stability analysis techniques are often employed to gain insight into the degree of stability of a system.

The most common methods used for determining the stability margins of a system, namely gain margin and phase margin, are based on frequency domain approaches such as the Nyquist, Bode, and Nichols method. These methods are limited to systems with less than two adjustable parameters. Systems that have more than two adjustable parameters require more complex control strategies such as the Vishnegradskii diagram, the parameter plane method or the stability equation method. These methods are used to plot the stability boundary of the system and to determine the effects of parameter variation on system stability but give no information about the stability margins (gain margins or phase margins) of the system.

The method described by Chang and Han [1] in their paper has been shown to be an effective method in determining the gain margin and phase margin of a system with adjustable parameters, such as a space shuttle. Their method combines commonly used frequency domain approaches with the parameter plane stability method to obtain boundary plots of constant gain margin and constant phase margin.

[1]  Chang, C-H., Han K-W. Gain Margins and Phase Margins for Control Systems with Adjustable Parameters. Journal of Guidance, Control, and Dynamics, (1989): 13(3), 404-408

$\mathrm{restart}\:$

The following section presents a brief overview of the approach taken by Chang and Han in determining the gain and phase margin of a system with adjustable parameters.

Consider the block diagram of a basic closed-loop control system shown in Figure 1.

Figure 1: Basic Block Diagram of a Closed-Loop System

The open loop transfer function for this system is defined as:

Substituting $s\=j\cdot \mathrm{\ω}$ into equation (1) yields:

Equation (2) can then be rewritten as:

Expressing equation (3) in terms of its magnitude and phase yields the following equation:

where

$∣G\left(\mathrm{j\ω}\right)∣\= \sqrt{{\left(\mathrm{\ℜG}\left(\mathrm{j\ω}\right)\right)}^2\+{\left(\mathrm{\ℑG}\left(\mathrm{j\ω}\right)\right)}^2}$

and

$\mathrm{\Φ}\mathit{ }\mathit{\=}\mathit{ }\angle G\left(\mathrm{j\ω}\right)\= \arctan \left\{\frac{\mathrm{\ℑG}\left(\mathrm{j\ω}\right)}{\mathrm{\ℜG}\left(\mathrm{j\ω}\right)}\right\}$

Combining equations (2) and (4) gives:

Dividing both sides of equation (5) by $∣G\left(\mathrm{j\ω}\right)∣{\ⅇ}^{j\cdot \mathrm{\Φ}}$ results in:  

Let's define:

$A\=\frac{1}{∣G\left(\mathrm{j\ω}\right)∣}$

and

$\mathrm{\Theta }\= \mathrm{\Phi } \+ 180$

where A is the gain margin of the system at $\mathrm{\Θ} \= 0$, and $\mathrm{\Theta }$ is the phase margin of the system when $A\=1$.

Then equation (6) can be rewritten as:

or as

The term $A{\ⅇ}^{-j\cdot \mathrm{\Θ}}$ is commonly referred to as the gain-phase margin tester. The gain-phase margin tester can be used to determine the gain and phase margin of a system with adjustable parameters by incorporating it as an additional block in the closed loop system, such as the one shown in Figure 2.

Figure 2: Basic Block Diagram of a Closed-Loop System with Gain-Phase Margin Tester

The gain-phase margin tester can also be expressed as: $$

where X and Y are:

Substituting X and Y into equation (8) yields:

which can be expressed in terms of the real and imaginary parts, namely:

Assuming that X and Y are parameters, we obtain the following expressions:

where $\mathrm{B\_\_1}$, $\mathrm{C\_\_1}$, $\mathrm{D\_\_1}$, $\mathrm{B\_\_2}$, $\mathrm{C\_\_2}$, and $\mathrm{D\_\_2}$ are functions of $\mathrm{\ω}$.

Solving equations (14) and (15) simultaneously yields:

and

where

If the system has adjustable parameters then equation (8) can be written as:

Assuming that $\mathrm{\alpha }$ and $\mathrm{\beta }$ are linear functions then

where $\mathrm{B\_\_1}$, $\mathrm{C\_\_1}$, $\mathrm{D\_\_1}$, $\mathrm{B\_\_2}$, $\mathrm{C\_\_2}$, and $\mathrm{D\_\_2}$ are functions of $\mathrm{\γ}\, ..\.\, A\, \mathrm{\Θ}\, \mathrm{\ω}$.

Solving equations (20) and (21) simultaneously yields:

where

If we let $A\=1$ and $\mathrm{\Θ}\=0$, and set $\mathrm{\gamma }$... equal to some constants, then as $\mathrm{\ω}$ varies from 0 to$\infty$ a locus that contains the stability boundary of the system without the gain-phase margin tester can be plotted in the $\mathrm{\alpha }$  vs. $\mathrm{\beta }$ plane. This locus can be considered the Nyquist plot of the system passing through the point $\left(1\,\mathrm{j0}\right)$.

If A is assumed equal to a constant value and $\mathrm{\Θ}\=0$, the locus formed in the $\mathrm{\alpha }$  vs. $\mathrm{\beta }$ plane is a boundary of constant gain margin. Similarly, if $\mathrm{\Θ}$ is assumed to be constant and $A\=1$, then the locus formed is a boundary of constant phase margin. The values of the phase crossover frequency and the gain crossover frequency can be determined by measuring the values of $\mathrm{\ω}$ on the gain margin boundary and the constant phase margin boundary.

The technique for determining gain and phase margin of a system with adjustable parameters, as elaborated upon in the previous section, will now be applied to the control system of an aerodynamic maneuverable re-entry vehicle.  The block diagram corresponding to the control system is shown in Figure 3.

Figure 3: Basic Block Diagram of an Aerodynamic Re-entry Vehicle

The transfer function associated with each block is defined below. The assumption is the system parameters, namely $\mathrm{\alpha }$, $\mathrm{\beta }$, $\mathrm{\Upsilon }$, and the frequency $\mathrm{\ω}$ at which the system is run are all real numbers.

$\mathrm{assume}\left(\mathrm{\alpha }\colon\colon \mathrm{real}\right) \:\mathrm{assume}\left(\mathrm{\beta }\colon\colon \mathrm{real}\right) \:\mathrm{assume}\left(\mathrm{\omega }\colon\colon \mathrm{real}\right) \:\mathrm{assume}\left(\mathrm{\Upsilon }\colon\colon \mathrm{real}\right)\:$

$\mathrm{G1}≔ \frac{Q}{\mathrm{x\_\_c}\'}\=\frac{7.7\cdot {\left(700\right)}^2}{s^2\+2\cdot \left(0.65\right)\cdot \left(700\right)\cdot s\+{\left(700\right)}^2}$

$\mathrm{G2}≔ \frac{\mathrm{x\_\_p}\'}{Q}\=\frac{\left(\frac{1}{2.2}\right)\cdot {\left(3030\right)}^2\cdot \left(s^2\+{5030}^2\right)}{{\left(5030\right)}^2\cdot \left(s^2\+{3030}^2\right)}$

$\mathrm{G3}≔\frac{\mathrm{x\_\_p}}{\mathrm{x\_\_p}\'}\=\frac{1}{s}$

$\mathrm{G4}≔\frac{\mathrm{x\_\_s}}{\mathrm{x\_\_p}}\=\frac{{5030}^2}{s^2\+{5030}^2}$

$\mathrm{G5}≔\frac{\mathrm{\θ}\mathrm{`\text{'}B`}}{\mathrm{x\_\_s}}\=\frac{\left(3.09\cdot {10}^{-4}\right)\cdot \left(s\+3.63\right)\cdot {\left(s^2\+{910}^2\right)}^{}}{s^2-{\left(24.8\right)}^2}$

$\mathrm{H1}≔\frac{\mathrm{y\_\_1}}{\mathrm{x\_\_p}}\=\frac{s}{s\+\mathrm{\ϒ}}$

$\mathrm{H2}≔ \mathrm{\α}$

$\mathrm{H3}≔\frac{\mathrm{\θ\_\_s}\'}{\mathrm{\θ\_\_B}\'}\=\frac{{\left(942\right)}^2}{s^2\+2\cdot 0.5\cdot 942\cdot s\+{942}^2}$

$\mathrm{H4}≔\mathrm{\β}$

The transfer function representing GH1 is:

$\mathrm{GH1}≔\mathrm{normal}\left(\frac{\left(\mathrm{rhs}\left(\mathrm{G1}\right)\cdot \mathrm{rhs}\left(\mathrm{G2}\right)\cdot \mathrm{rhs}\left(\mathrm{G3}\right)\right)}{1\+\left(\mathrm{rhs}\left(\mathrm{G1}\right)\cdot \mathrm{rhs}\left(\mathrm{G2}\right)\cdot \mathrm{rhs}\left(\mathrm{G3}\right)\right)\cdot \left(\mathrm{rhs}\left(\mathrm{H1}\right)\cdot \mathrm{H2}\right)}\right)$

Similarly, the transfer function for GH2 is:

$\mathrm{GH2}≔\mathrm{rhs}\left(\mathrm{G4}\right)\cdot \mathrm{rhs}\left(\mathrm{G5}\right)\cdot \mathrm{rhs}\left(\mathrm{H3}\right)\cdot \mathrm{H4}$

The transfer function for the simplified model is found by combining the transfer function equations for GH1 and GH2.

$\mathrm{GHsimplified}≔\mathrm{GH1}\cdot \mathrm{GH2}$

To determine the stability of the system, the closed-loop system equation as defined in $\mathrm{GHsimplified}$ was modified to accommodate the inclusion of the gain-phase margin tester as defined in equation (9). This results in an equation of the form:

$G\=\mathrm{Denominator}\+\left(A\cdot \cos \left(\mathrm{\Theta }\right)-j\cdot A\cdot \sin \left(\mathrm{\Theta }\right)\right)\cdot \mathrm{Numerator}\:$

where Denominator and Numerator refer to the denominator and numerator of $\mathrm{GHsimplified}$, respectively. The numerator and denominator of $\mathrm{GHsimplified}$ can be extracted using the $\mathrm{numer}$ and $\mathrm{denom}$ commands.

$\mathrm{NumerGHsimplified}≔\mathrm{expand}\left(\mathrm{numer}\left(\mathrm{GHsimplified}\right)\,s\right)$

$\mathrm{DenomGHsimplified}≔\mathrm{expand}\left(\mathrm{denom}\left(\mathrm{GHsimplified}\right)\,s\right)$

Using equations (37) and (38) the closed-loop system equation with the addition of a gain-phase margin tester block can now be obtained.

$\mathrm{GH\_\_GainPhaseMarginTester}≔\mathrm{DenomGHsimplified}\+\left(A\cdot \cos \left(\mathrm{\Theta }\right)-j\cdot A\cdot \sin \left(\mathrm{\Theta }\right)\right)\cdot \mathrm{NumerGHsimplified}$

The frequency response of the modified closed-loop system is obtained by replacing the term s with $\mathrm{j\ω}$.

$\mathrm{GH\ω\_\_GainPhaseMarginTester}≔\mathrm{expand}\left(\mathrm{subs}\left(s\=j\cdot \mathrm{\ω}\,\mathrm{GH\_\_GainPhaseMarginTester}\right)\right)$

To see the effects of parameter variations on the system, the equations ,  and  defined respectively in (22), (23) and (24), will be calculated.

Recalling equations (20) and (21), respectively  and , note that:

The real coefficients corresponding to  are extracted using the following commands:

$\mathrm{RealValues}≔\mathrm{seq}\left(\mathrm{map}\left(\mathrm{evalc}\,\mathrm{\Re }\left(\mathrm{op}\left(\mathrm{GH\ω\_\_GainPhaseMarginTester}\right)\left[\mathrm{JJ}\right]\right)\right)\,\mathrm{JJ}\=1..\mathrm{nops}\left(\mathrm{GH\ω\_\_GainPhaseMarginTester}\right)\right)\:$ 

$\mathrm{RealCoeff}≔\mathrm{coeffs}\left(\mathrm{add}\left(\mathrm{RealValues}\left[\mathrm{JJ}\right]\,\mathrm{JJ}\=1..\mathrm{nops}\left(\mathrm{GH\ω\_\_GainPhaseMarginTester}\right)\right)\,\left[\mathrm{\α}\,\mathrm{\β}\right]\right)\:$

Similarly, the imaginary coefficients corresponding to  are extracted using the following commands:

$\mathrm{ImaginaryValues}≔\mathrm{seq}\left(\mathrm{map}\left(\mathrm{evalc}\,\mathrm{\Im }\left(\mathrm{op}\left(\mathrm{GH\ω\_\_GainPhaseMarginTester}\right)\left[\mathrm{JJ}\right]\right)\right)\,\mathrm{JJ}\=1..\mathrm{nops}\left(\mathrm{GH\ω\_\_GainPhaseMarginTester}\right)\right)\:$

$\mathrm{ImaginaryCoeff}≔\mathrm{coeffs}\left(\mathrm{add}\left(\mathrm{ImaginaryValues}\left[\mathrm{JJ}\right]\,\mathrm{JJ}\=1..\mathrm{nops}\left(\mathrm{GH\ω\_\_GainPhaseMarginTester}\right)\right)\,\left[\mathrm{\α}\,\mathrm{\β}\right]\right)\:$

Using the sort command the values for $\mathrm{B1}$, $\mathrm{C1}$, $\mathrm{D1}$, $\mathrm{B2}$, $\mathrm{C2}$, $\mathrm{D2}$, $\mathrm{\Δ1}$ and $\mathrm{\Δ2}$ can be easily obtained.

$\mathrm{B1}≔\mathrm{sort}\left(\mathrm{RealCoeff}\left[2\right]\,\mathrm{\ω}\right)$

$\mathrm{C1}≔\mathrm{sort}\left(\mathrm{RealCoeff}\left[3\right]\,\mathrm{\ω}\right)$

$\mathrm{D1}≔\mathrm{sort}\left(\mathrm{RealCoeff}\left[1\right]\,\mathrm{\ω}\right)$

$\mathrm{B2}≔\mathrm{sort}\left(\mathrm{ImaginaryCoeff}\left[2\right]\,\mathrm{\ω}\right)$

$\mathrm{C2}≔\mathrm{sort}\left(\mathrm{ImaginaryCoeff}\left[3\right]\,\mathrm{\ω}\right)$

$\mathrm{D2}≔\mathrm{sort}\left(\mathrm{ImaginaryCoeff}\left[1\right]\,\mathrm{\ω}\right)$

$\mathrm{\Δ}≔\mathrm{B1}\cdot \mathrm{C2}-\mathrm{B2}\cdot \mathrm{C1}$

Using the values obtained for $\mathrm{B1}\, \mathrm{C1}\, \mathrm{D1}\, \mathrm{B2}\, \mathrm{C2}\, \mathrm{D2}\, \mathrm{\Δ1}$ and $\mathrm{\Δ2}$, the stability boundary plots for the system in the $\mathrm{\alpha }$ vs. $\mathrm{\beta }$ plane can be obtained for different values of $\mathrm{\α}$ and $\mathrm{\β}$. For the sake of efficiency, Maple procedures to calculate the values of $\mathrm{\α}$ and $\mathrm{\β}$are generated from the equations for $\mathrm{\α}$ and $\mathrm{\β}$ as defined in equations (22) and (23).

$\mathrm{\α\_proc}≔\mathrm{codegen}\left[\mathrm{makeproc}\right]\left(\mathrm{rhs}\left(\right)\,\mathrm{parameters}\=\left[\mathrm{\ϒ}\colon\colon \mathrm{numeric}\,A\mathit{\colon\colon }\mathrm{numeric}\,\mathrm{\Theta }\colon\colon \mathrm{numeric}\,\mathrm{\omega }\colon\colon \mathrm{numeric}\right]\right)\:$

$\mathrm{\α\_optimize}≔\mathrm{codegen}\left[\mathrm{optimize}\right]\left(\mathrm{\α\_proc}\right)\:$

$\mathrm{\β\_proc}≔\mathrm{codegen}\left[\mathrm{makeproc}\right]\left(\mathrm{rhs}\left(\right)\mathit{\,}\mathrm{parameters}\=\left[\mathrm{\Upsilon }\mathit{\colon\colon }\mathrm{numeric}\,A\mathit{\colon\colon }\mathrm{numeric}\mathit{\,}\mathrm{\Θ}\mathit{\colon\colon }\mathrm{numeric}\mathit{\,}\mathrm{\ω}\mathit{\colon\colon }\mathrm{numeric}\right]\right)\:$

$\mathrm{\β\_optimize}≔\mathrm{codegen}\left[\mathrm{optimize}\right]\left(\mathrm{\β\_proc}\right)\:$

The stability curves shown below assume the value of $\mathrm{\ϒ}$ to be 30. The Maple code used to generated the plots is contained in the following code edit region.

Combining the boundary plots of constant phase margins and constant gain margins that lie within the stable region yields the plot shown below. Each region has a specified gain and phase margin. For instance, the region defined by $\mathrm{M\_\_1}$ will have a gain margin of $A\&gt;3$ and $A<1\&sol;3$, and a phase margin of $30°<\mathrm{\&Theta;} <45°$. While the region as defined by $\mathrm{M\_\_2}$ will have a gain margin of $A\&gt;3$ and $1\&sol;3<A<1\&sol;2$, and a phase margin of $15°<\mathrm{\&Theta;} <30°$.

The phase crossover frequency values can be obtained for any point along the constant gain margin boundary curves. Similarly, the gain crossover frequency values can be obtained for any point along the constant phase margin boundary curves. The phase crossover frequency values and the gain crossover frequency values for 6 data points are listed in the table below. The code used to obtain these values can be found in the following code edit region.