# Inner/outer control design for vertical takeoff and landing aircraft¶

This script demonstrates the use of the python-control package for analysis and design of a controller for a vectored thrust aircraft model that is used as a running example through the text Feedback Systems by Astrom and Murray. This example makes use of MATLAB compatible commands.

## Code¶

```  1# pvtol-nested.py - inner/outer design for vectored thrust aircraft
2# RMM, 5 Sep 09
3#
4# This file works through a fairly complicated control design and
5# analysis, corresponding to the planar vertical takeoff and landing
6# (PVTOL) aircraft in Astrom and Murray, Chapter 11.  It is intended
7# to demonstrate the basic functionality of the python-control
8# package.
9#
10
11import os
12import matplotlib.pyplot as plt  # MATLAB plotting functions
13from control.matlab import *    # MATLAB-like functions
14import numpy as np
15
16# System parameters
17m = 4               # mass of aircraft
18J = 0.0475          # inertia around pitch axis
19r = 0.25            # distance to center of force
20g = 9.8             # gravitational constant
21c = 0.05            # damping factor (estimated)
22
23# Transfer functions for dynamics
24Pi = tf([r], [J, 0, 0])  # inner loop (roll)
25Po = tf([1], [m, c, 0])  # outer loop (position)
26
27#
28# Inner loop control design
29#
30# This is the controller for the pitch dynamics.  Goal is to have
31# fast response for the pitch dynamics so that we can use this as a
32# control for the lateral dynamics
33#
34
35# Design a simple lead controller for the system
36k, a, b = 200, 2, 50
37Ci = k*tf([1, a], [1, b])  # lead compensator
38Li = Pi*Ci
39
40# Bode plot for the open loop process
41plt.figure(1)
42bode(Pi)
43
44# Bode plot for the loop transfer function, with margins
45plt.figure(2)
46bode(Li)
47
48# Compute out the gain and phase margins
49#! Not implemented
50# gm, pm, wcg, wcp = margin(Li)
51
52# Compute the sensitivity and complementary sensitivity functions
53Si = feedback(1, Li)
54Ti = Li*Si
55
56# Check to make sure that the specification is met
57plt.figure(3)
58gangof4(Pi, Ci)
59
60# Compute out the actual transfer function from u1 to v1 (see L8.2 notes)
61# Hi = Ci*(1-m*g*Pi)/(1+Ci*Pi)
62Hi = parallel(feedback(Ci, Pi), -m*g*feedback(Ci*Pi, 1))
63
64plt.figure(4)
65plt.clf()
66plt.subplot(221)
67bode(Hi)
68
69# Now design the lateral control system
70a, b, K = 0.02, 5, 2
71Co = -K*tf([1, 0.3], [1, 10])  # another lead compensator
72Lo = -m*g*Po*Co
73
74plt.figure(5)
75bode(Lo)  # margin(Lo)
76
77# Finally compute the real outer-loop loop gain + responses
78L = Co*Hi*Po
79S = feedback(1, L)
80T = feedback(L, 1)
81
82# Compute stability margins
83gm, pm, wgc, wpc = margin(L)
84print("Gain margin: %g at %g" % (gm, wgc))
85print("Phase margin: %g at %g" % (pm, wpc))
86
87plt.figure(6)
88plt.clf()
89bode(L, np.logspace(-4, 3))
90
91# Add crossover line to the magnitude plot
92#
93# Note: in matplotlib before v2.1, the following code worked:
94#
95#   plt.subplot(211); hold(True);
96#   loglog([1e-4, 1e3], [1, 1], 'k-')
97#
98# In later versions of matplotlib the call to plt.subplot will clear the
99# axes and so we have to extract the axes that we want to use by hand.
100# In addition, hold() is deprecated so we no longer require it.
101#
102for ax in plt.gcf().axes:
103    if ax.get_label() == 'control-bode-magnitude':
104        break
105ax.semilogx([1e-4, 1e3], 20*np.log10([1, 1]), 'k-')
106
107#
108# Replot phase starting at -90 degrees
109#
110# Get the phase plot axes
111for ax in plt.gcf().axes:
112    if ax.get_label() == 'control-bode-phase':
113        break
114
115# Recreate the frequency response and shift the phase
116mag, phase, w = freqresp(L, np.logspace(-4, 3))
117phase = phase - 360
118
119# Replot the phase by hand
120ax.semilogx([1e-4, 1e3], [-180, -180], 'k-')
121ax.semilogx(w, np.squeeze(phase), 'b-')
122ax.axis([1e-4, 1e3, -360, 0])
123plt.xlabel('Frequency [deg]')
124plt.ylabel('Phase [deg]')
125# plt.set(gca, 'YTick', [-360, -270, -180, -90, 0])
126# plt.set(gca, 'XTick', [10^-4, 10^-2, 1, 100])
127
128#
129# Nyquist plot for complete design
130#
131plt.figure(7)
132plt.clf()
133nyquist(L, (0.0001, 1000))
134
135# Add a box in the region we are going to expand
136plt.plot([-2, -2, 1, 1, -2], [-4, 4, 4, -4, -4], 'r-')
137
138# Expanded region
139plt.figure(8)
140plt.clf()
141nyquist(L)
142plt.axis([-2, 1, -4, 4])
143
144# set up the color
145color = 'b'
146
147# Add arrows to the plot
148# H1 = L.evalfr(0.4); H2 = L.evalfr(0.41);
149# arrow([real(H1), imag(H1)], [real(H2), imag(H2)], AM_normal_arrowsize, \
150#  'EdgeColor', color, 'FaceColor', color);
151
152# H1 = freqresp(L, 0.35); H2 = freqresp(L, 0.36);
153# arrow([real(H2), -imag(H2)], [real(H1), -imag(H1)], AM_normal_arrowsize, \
154#  'EdgeColor', color, 'FaceColor', color);
155
156plt.figure(9)
157Yvec, Tvec = step(T, np.linspace(0, 20))
158plt.plot(Tvec.T, Yvec.T)
159
160Yvec, Tvec = step(Co*S, np.linspace(0, 20))
161plt.plot(Tvec.T, Yvec.T)
162
163plt.figure(10)
164plt.clf()
165P, Z = pzmap(T, plot=True, grid=True)
166print("Closed loop poles and zeros: ", P, Z)
167
168# Gang of Four
169plt.figure(11)
170plt.clf()
171gangof4(Hi*Po, Co)
172
173if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
174    plt.show()
```

## Notes¶

1. The environment variable PYCONTROL_TEST_EXAMPLES is used for testing to turn off plotting of the outputs.