# Model-Reference Adaptive Control (MRAC) SISO, direct Lyapunov rule

## Code

  1# mrac_siso_lyapunov.py
2# Johannes Kaisinger, 3 July 2023
3#
4# Demonstrate a MRAC example for a SISO plant using Lyapunov rule.
5# Based on [1] Ex 5.7, Fig 5.12 & 5.13.
6# Notation as in [2].
7#
8# [1] K. J. Aström & B. Wittenmark "Adaptive Control" Second Edition, 2008.
9#
10# [2] Nhan T. Nguyen "Model-Reference Adaptive Control", 2018.
11
12import numpy as np
13import scipy.signal as signal
14import matplotlib.pyplot as plt
15import os
16
17import control as ct
18
19# Plant model as linear state-space system
20A = -1
21B = 0.5
22C = 1
23D = 0
24
25io_plant = ct.ss(A, B, C, D,
26                 inputs=('u'), outputs=('x'), states=('x'), name='plant')
27
28# Reference model as linear state-space system
29Am = -2
30Bm = 2
31Cm = 1
32Dm = 0
33
34io_ref_model = ct.ss(Am, Bm, Cm, Dm,
35                     inputs=('r'), outputs=('xm'), states=('xm'), name='ref_model')
36
37# Adaptive control law, u = kx*x + kr*r
38kr_star = (Bm)/B
39print(f"Optimal value for {kr_star = }")
40kx_star = (Am-A)/B
41print(f"Optimal value for {kx_star = }")
42
44    """Internal state of adaptive controller, f(t,x,u;p)"""
45
46    # Parameters
47    gam = params["gam"]
48    signB = params["signB"]
49
50    # Controller inputs
51    r = uc[0]
52    xm = uc[1]
53    x = uc[2]
54
55    # Controller states
56    # x1 = xc[0] # kr
57    # x2 = xc[1] # kx
58
59    # Algebraic relationships
60    e = xm - x
61
62    # Controller dynamics
63    d_x1 = gam*r*e*signB
64    d_x2 = gam*x*e*signB
65
66    return [d_x1, d_x2]
67
69    """Algebraic output from adaptive controller, g(t,x,u;p)"""
70
71    # Controller inputs
72    r = uc[0]
73    #xm = uc[1]
74    x = uc[2]
75
76    # Controller state
77    kr = xc[0]
78    kx = xc[1]
79
80    # Control law
81    u = kx*x + kr*r
82
83    return [u]
84
85params={"gam":1, "Am":Am, "Bm":Bm, "signB":np.sign(B)}
86
87io_controller = ct.nlsys(
90    inputs=('r', 'xm', 'x'),
91    outputs=('u'),
92    states=2,
93    params=params,
94    name='control',
95    dt=0
96)
97
98# Overall closed loop system
99io_closed = ct.interconnect(
100    [io_plant, io_ref_model, io_controller],
101    connections=[
102        ['plant.u', 'control.u'],
103        ['control.xm', 'ref_model.xm'],
104        ['control.x', 'plant.x']
105    ],
106    inplist=['control.r', 'ref_model.r'],
107    outlist=['plant.x', 'control.u'],
108    dt=0
109)
110
111# Set simulation duration and time steps
112Tend = 100
113dt = 0.1
114
115# Define simulation time
116t_vec = np.arange(0, Tend, dt)
117
118# Define control reference input
119r_vec = np.zeros((2, len(t_vec)))
120rect = signal.square(2 * np.pi * 0.05 * t_vec)
121r_vec[0, :] = rect
122r_vec[1, :] = r_vec[0, :]
123
124plt.figure(figsize=(16,8))
125plt.plot(t_vec, r_vec[0,:])
126plt.title(r'reference input $r$')
127plt.show()
128
129# Set initial conditions, io_closed
130X0 = np.zeros((4, 1))
131X0[0] = 0 # state of plant, (x)
132X0[1] = 0 # state of ref_model, (xm)
133X0[2] = 0 # state of controller, (kr)
134X0[3] = 0 # state of controller, (kx)
135
136# Simulate the system with different gammas
137tout1, yout1, xout1 = ct.input_output_response(io_closed, t_vec, r_vec, X0,
138                                               return_x=True, params={"gam":0.2})
139tout2, yout2, xout2 = ct.input_output_response(io_closed, t_vec, r_vec, X0,
140                                               return_x=True, params={"gam":1.0})
141tout3, yout3, xout3 = ct.input_output_response(io_closed, t_vec, r_vec, X0,
142                                               return_x=True, params={"gam":5.0})
143
144plt.figure(figsize=(16,8))
145plt.subplot(2,1,1)
146plt.plot(tout1, yout1[0,:], label=r'$x_{\gamma = 0.2}$')
147plt.plot(tout2, yout2[0,:], label=r'$x_{\gamma = 1.0}$')
148plt.plot(tout2, yout3[0,:], label=r'$x_{\gamma = 5.0}$')
149plt.plot(tout1, xout1[1,:] ,label=r'$x_{m}$', color='black', linestyle='--')
150plt.legend(fontsize=14)
151plt.title(r'system response $x, (x_m)$')
152plt.subplot(2,1,2)
153plt.plot(tout1, yout1[1,:], label=r'$u_{\gamma = 0.2}$')
154plt.plot(tout2, yout2[1,:], label=r'$u_{\gamma = 1.0}$')
155plt.plot(tout3, yout3[1,:], label=r'$u_{\gamma = 5.0}$')
156plt.legend(loc=4, fontsize=14)
157plt.title(r'control $u$')
158
159plt.figure(figsize=(16,8))
160plt.subplot(2,1,1)
161plt.plot(tout1, xout1[2,:], label=r'$k_{r, \gamma = 0.2}$')
162plt.plot(tout2, xout2[2,:], label=r'$k_{r, \gamma = 1.0}$')
163plt.plot(tout3, xout3[2,:], label=r'$k_{r, \gamma = 5.0}$')
164plt.hlines(kr_star, 0, Tend, label=r'$k_r^{\ast}$', color='black', linestyle='--')
165plt.legend(loc=4, fontsize=14)
166plt.title(r'control gain $k_r$ (feedforward)')
167plt.subplot(2,1,2)
168plt.plot(tout1, xout1[3,:], label=r'$k_{x, \gamma = 0.2}$')
169plt.plot(tout2, xout2[3,:], label=r'$k_{x, \gamma = 1.0}$')
170plt.plot(tout3, xout3[3,:], label=r'$k_{x, \gamma = 5.0}$')
171plt.hlines(kx_star, 0, Tend, label=r'$k_x^{\ast}$', color='black', linestyle='--')
172plt.legend(loc=4, fontsize=14)
173plt.title(r'control gain $k_x$ (feedback)')
174if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
175    plt.show()


## Notes

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