Model-Reference Adaptive Control (MRAC) SISO, direct MIT rule

Code

  1# mrac_siso_mit.py
  2# Johannes Kaisinger, 3 July 2023
  3#
  4# Demonstrate a MRAC example for a SISO plant using MIT rule.
  5# Based on [1] Ex 5.2, Fig 5.5 & 5.6.
  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
 43def adaptive_controller_state(t, xc, uc, params):
 44    """Internal state of adaptive controller, f(t,x,u;p)"""
 45    
 46    # Parameters
 47    gam = params["gam"]
 48    Am = params["Am"]
 49    Bm = params["Bm"]
 50    signB = params["signB"]
 51
 52    # Controller inputs
 53    r = uc[0]
 54    xm = uc[1]
 55    x = uc[2]
 56
 57    # Controller states
 58    x1 = xc[0] #
 59    # x2 = xc[1] # kr
 60    x3 = xc[2] #
 61    # x4 = xc[3] # kx
 62    
 63    # Algebraic relationships
 64    e = xm - x
 65
 66    # Controller dynamics
 67    d_x1 = Am*x1 + Am*r
 68    d_x2 = - gam*x1*e*signB
 69    d_x3 = Am*x3 + Am*x
 70    d_x4 = - gam*x3*e*signB
 71
 72    return [d_x1, d_x2, d_x3, d_x4]
 73
 74def adaptive_controller_output(t, xc, uc, params):
 75    """Algebraic output from adaptive controller, g(t,x,u;p)"""
 76
 77    # Controller inputs
 78    r = uc[0]
 79    # xm = uc[1]
 80    x = uc[2]
 81    
 82    # Controller state
 83    kr = xc[1]
 84    kx = xc[3]
 85    
 86    # Control law
 87    u = kx*x + kr*r
 88
 89    return [u]
 90
 91params={"gam":1, "Am":Am, "Bm":Bm, "signB":np.sign(B)}
 92
 93io_controller = ct.nlsys(
 94    adaptive_controller_state,
 95    adaptive_controller_output,
 96    inputs=('r', 'xm', 'x'),
 97    outputs=('u'),
 98    states=4,
 99    params=params,
100    name='control',
101    dt=0
102)
103
104# Overall closed loop system
105io_closed = ct.interconnect(
106    [io_plant, io_ref_model, io_controller],
107    connections=[
108        ['plant.u', 'control.u'],
109        ['control.xm', 'ref_model.xm'],
110        ['control.x', 'plant.x']
111    ],
112    inplist=['control.r', 'ref_model.r'],
113    outlist=['plant.x', 'control.u'],
114    dt=0
115)
116
117# Set simulation duration and time steps
118Tend = 100
119dt = 0.1
120
121# Define simulation time 
122t_vec = np.arange(0, Tend, dt)
123
124# Define control reference input
125r_vec = np.zeros((2, len(t_vec)))
126square = signal.square(2 * np.pi * 0.05 * t_vec)
127r_vec[0, :] = square
128r_vec[1, :] = r_vec[0, :]
129
130plt.figure(figsize=(16,8))
131plt.plot(t_vec, r_vec[0,:])
132plt.title(r'reference input $r$')
133plt.show()
134
135# Set initial conditions, io_closed
136X0 = np.zeros((6, 1))
137X0[0] = 0 # state of plant, (x)
138X0[1] = 0 # state of ref_model, (xm)
139X0[2] = 0 # state of controller,
140X0[3] = 0 # state of controller, (kr)
141X0[4] = 0 # state of controller,
142X0[5] = 0 # state of controller, (kx)
143
144# Simulate the system with different gammas
145tout1, yout1, xout1 = ct.input_output_response(io_closed, t_vec, r_vec, X0,
146                                               return_x=True, params={"gam":0.2})
147tout2, yout2, xout2 = ct.input_output_response(io_closed, t_vec, r_vec, X0,
148                                               return_x=True, params={"gam":1.0})
149tout3, yout3, xout3 = ct.input_output_response(io_closed, t_vec, r_vec, X0,
150                                               return_x=True, params={"gam":5.0})
151
152plt.figure(figsize=(16,8))
153plt.subplot(2,1,1)
154plt.plot(tout1, yout1[0,:], label=r'$x_{\gamma = 0.2}$')
155plt.plot(tout2, yout2[0,:], label=r'$x_{\gamma = 1.0}$')
156plt.plot(tout2, yout3[0,:], label=r'$x_{\gamma = 5.0}$')
157plt.plot(tout1, xout1[1,:] ,label=r'$x_{m}$', color='black', linestyle='--')
158plt.legend(fontsize=14)
159plt.title(r'system response $x, (x_m)$')
160plt.subplot(2,1,2)
161plt.plot(tout1, yout1[1,:], label=r'$u_{\gamma = 0.2}$')
162plt.plot(tout2, yout2[1,:], label=r'$u_{\gamma = 1.0}$')
163plt.plot(tout3, yout3[1,:], label=r'$u_{\gamma = 5.0}$')
164plt.legend(loc=4, fontsize=14)
165plt.title(r'control $u$')
166
167plt.figure(figsize=(16,8))
168plt.subplot(2,1,1)
169plt.plot(tout1, xout1[3,:], label=r'$k_{r, \gamma = 0.2}$')
170plt.plot(tout2, xout2[3,:], label=r'$k_{r, \gamma = 1.0}$')
171plt.plot(tout3, xout3[3,:], label=r'$k_{r, \gamma = 5.0}$')
172plt.hlines(kr_star, 0, Tend, label=r'$k_r^{\ast}$', color='black', linestyle='--')
173plt.legend(loc=4, fontsize=14)
174plt.title(r'control gain $k_r$ (feedforward)')
175plt.subplot(2,1,2)
176plt.plot(tout1, xout1[5,:], label=r'$k_{x, \gamma = 0.2}$')
177plt.plot(tout2, xout2[5,:], label=r'$k_{x, \gamma = 1.0}$')
178plt.plot(tout3, xout3[5,:], label=r'$k_{x, \gamma = 5.0}$')
179plt.hlines(kx_star, 0, Tend, label=r'$k_x^{\ast}$', color='black', linestyle='--')
180plt.legend(loc=4, fontsize=14)
181plt.title(r'control gain $k_x$ (feedback)')
182
183if 'PYCONTROL_TEST_EXAMPLES' not in os.environ:
184    plt.show()

Notes

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