From f3713f1f9bd5e399a28b88ae773654538b5cc65a Mon Sep 17 00:00:00 2001 From: Johannes Kaisinger Date: Fri, 16 Jun 2023 00:43:13 +0200 Subject: [PATCH] Add direct mrac siso examples, using new input/output system (ss, nlsys, interconnect) --- doc/examples.rst | 2 + doc/mrac_siso_lyapunov.py | 1 + doc/mrac_siso_lyapunov.rst | 15 +++ doc/mrac_siso_mit.py | 1 + doc/mrac_siso_mit.rst | 15 +++ examples/mrac_siso_lyapunov.py | 174 +++++++++++++++++++++++++++++++ examples/mrac_siso_mit.py | 182 +++++++++++++++++++++++++++++++++ 7 files changed, 390 insertions(+) create mode 120000 doc/mrac_siso_lyapunov.py create mode 100644 doc/mrac_siso_lyapunov.rst create mode 120000 doc/mrac_siso_mit.py create mode 100644 doc/mrac_siso_mit.rst create mode 100644 examples/mrac_siso_lyapunov.py create mode 100644 examples/mrac_siso_mit.py diff --git a/doc/examples.rst b/doc/examples.rst index 1d8ac3f42..41e2b42d6 100644 --- a/doc/examples.rst +++ b/doc/examples.rst @@ -33,6 +33,8 @@ other sources. steering-gainsched steering-optimal kincar-flatsys + mrac_siso_mit + mrac_siso_lyapunov Jupyter notebooks ================= diff --git a/doc/mrac_siso_lyapunov.py b/doc/mrac_siso_lyapunov.py new file mode 120000 index 000000000..aaccf5585 --- /dev/null +++ b/doc/mrac_siso_lyapunov.py @@ -0,0 +1 @@ +../examples/mrac_siso_lyapunov.py \ No newline at end of file diff --git a/doc/mrac_siso_lyapunov.rst b/doc/mrac_siso_lyapunov.rst new file mode 100644 index 000000000..525968882 --- /dev/null +++ b/doc/mrac_siso_lyapunov.rst @@ -0,0 +1,15 @@ +Model-Reference Adaptive Control (MRAC) SISO, direct Lyapunov rule +------------------------------------------------------------------ + +Code +.... +.. literalinclude:: mrac_siso_lyapunov.py + :language: python + :linenos: + + +Notes +..... + +1. The environment variable `PYCONTROL_TEST_EXAMPLES` is used for +testing to turn off plotting of the outputs. \ No newline at end of file diff --git a/doc/mrac_siso_mit.py b/doc/mrac_siso_mit.py new file mode 120000 index 000000000..b6a226f7c --- /dev/null +++ b/doc/mrac_siso_mit.py @@ -0,0 +1 @@ +../examples/mrac_siso_mit.py \ No newline at end of file diff --git a/doc/mrac_siso_mit.rst b/doc/mrac_siso_mit.rst new file mode 100644 index 000000000..8be834d6d --- /dev/null +++ b/doc/mrac_siso_mit.rst @@ -0,0 +1,15 @@ +Model-Reference Adaptive Control (MRAC) SISO, direct MIT rule +------------------------------------------------------------- + +Code +.... +.. literalinclude:: mrac_siso_mit.py + :language: python + :linenos: + + +Notes +..... + +1. The environment variable `PYCONTROL_TEST_EXAMPLES` is used for +testing to turn off plotting of the outputs.0 \ No newline at end of file diff --git a/examples/mrac_siso_lyapunov.py b/examples/mrac_siso_lyapunov.py new file mode 100644 index 000000000..00dbf63aa --- /dev/null +++ b/examples/mrac_siso_lyapunov.py @@ -0,0 +1,174 @@ +# mrac_siso_lyapunov.py +# Johannes Kaisinger, 3 July 2023 +# +# Demonstrate a MRAC example for a SISO plant using Lyapunov rule. +# Based on [1] Ex 5.7, Fig 5.12 & 5.13. +# Notation as in [2]. +# +# [1] K. J. Aström & B. Wittenmark "Adaptive Control" Second Edition, 2008. +# +# [2] Nhan T. Nguyen "Model-Reference Adaptive Control", 2018. + +import numpy as np +import scipy.signal as signal +import matplotlib.pyplot as plt + +import control as ct + +# Plant model as linear state-space system +A = -1 +B = 0.5 +C = 1 +D = 0 + +io_plant = ct.ss(A, B, C, D, + inputs=('u'), outputs=('x'), states=('x'), name='plant') + +# Reference model as linear state-space system +Am = -2 +Bm = 2 +Cm = 1 +Dm = 0 + +io_ref_model = ct.ss(Am, Bm, Cm, Dm, + inputs=('r'), outputs=('xm'), states=('xm'), name='ref_model') + +# Adaptive control law, u = kx*x + kr*r +kr_star = (Bm)/B +print(f"Optimal value for {kr_star = }") +kx_star = (Am-A)/B +print(f"Optimal value for {kx_star = }") + +def adaptive_controller_state(_t, xc, uc, params): + """Internal state of adaptive controller, f(t,x,u;p)""" + + # Parameters + gam = params["gam"] + signB = params["signB"] + + # Controller inputs + r = uc[0] + xm = uc[1] + x = uc[2] + + # Controller states + # x1 = xc[0] # kr + # x2 = xc[1] # kx + + # Algebraic relationships + e = xm - x + + # Controller dynamics + d_x1 = gam*r*e*signB + d_x2 = gam*x*e*signB + + return [d_x1, d_x2] + +def adaptive_controller_output(_t, xc, uc, params): + """Algebraic output from adaptive controller, g(t,x,u;p)""" + + # Controller inputs + r = uc[0] + #xm = uc[1] + x = uc[2] + + # Controller state + kr = xc[0] + kx = xc[1] + + # Control law + u = kx*x + kr*r + + return [u] + +params={"gam":1, "Am":Am, "Bm":Bm, "signB":np.sign(B)} + +io_controller = ct.nlsys( + adaptive_controller_state, + adaptive_controller_output, + inputs=('r', 'xm', 'x'), + outputs=('u'), + states=2, + params=params, + name='control', + dt=0 +) + +# Overall closed loop system +io_closed = ct.interconnect( + [io_plant, io_ref_model, io_controller], + connections=[ + ['plant.u', 'control.u'], + ['control.xm', 'ref_model.xm'], + ['control.x', 'plant.x'] + ], + inplist=['control.r', 'ref_model.r'], + outlist=['plant.x', 'control.u'], + dt=0 +) + +# Set simulation duration and time steps +Tend = 100 +dt = 0.1 + +# Define simulation time +t_vec = np.arange(0, Tend, dt) + +# Define control reference input +r_vec = np.zeros((2, len(t_vec))) +rect = signal.square(2 * np.pi * 0.05 * t_vec) +r_vec[0, :] = rect +r_vec[1, :] = r_vec[0, :] + +plt.figure(figsize=(16,8)) +plt.plot(t_vec, r_vec[0,:]) +plt.title(r'reference input $r$') +plt.show() + +# Set initial conditions, io_closed +X0 = np.zeros((4, 1)) +X0[0] = 0 # state of plant, (x) +X0[1] = 0 # state of ref_model, (xm) +X0[2] = 0 # state of controller, (kr) +X0[3] = 0 # state of controller, (kx) + +# Simulate the system with different gammas +tout1, yout1, xout1 = ct.input_output_response(io_closed, t_vec, r_vec, X0, + return_x=True, params={"gam":0.2}) +tout2, yout2, xout2 = ct.input_output_response(io_closed, t_vec, r_vec, X0, + return_x=True, params={"gam":1.0}) +tout3, yout3, xout3 = ct.input_output_response(io_closed, t_vec, r_vec, X0, + return_x=True, params={"gam":5.0}) + +plt.figure(figsize=(16,8)) +plt.subplot(2,1,1) +plt.plot(tout1, yout1[0,:], label=r'$x_{\gamma = 0.2}$') +plt.plot(tout2, yout2[0,:], label=r'$x_{\gamma = 1.0}$') +plt.plot(tout2, yout3[0,:], label=r'$x_{\gamma = 5.0}$') +plt.plot(tout1, xout1[1,:] ,label=r'$x_{m}$', color='black', linestyle='--') +plt.legend(fontsize=14) +plt.title(r'system response $x, (x_m)$') +plt.subplot(2,1,2) +plt.plot(tout1, yout1[1,:], label=r'$u_{\gamma = 0.2}$') +plt.plot(tout2, yout2[1,:], label=r'$u_{\gamma = 1.0}$') +plt.plot(tout3, yout3[1,:], label=r'$u_{\gamma = 5.0}$') +plt.legend(loc=4, fontsize=14) +plt.title(r'control $u$') +plt.show() + +plt.figure(figsize=(16,8)) +plt.subplot(2,1,1) +plt.plot(tout1, xout1[2,:], label=r'$k_{r, \gamma = 0.2}$') +plt.plot(tout2, xout2[2,:], label=r'$k_{r, \gamma = 1.0}$') +plt.plot(tout3, xout3[2,:], label=r'$k_{r, \gamma = 5.0}$') +plt.hlines(kr_star, 0, Tend, label=r'$k_r^{\ast}$', color='black', linestyle='--') +plt.legend(loc=4, fontsize=14) +plt.title(r'control gain $k_r$ (feedforward)') +plt.subplot(2,1,2) +plt.plot(tout1, xout1[3,:], label=r'$k_{x, \gamma = 0.2}$') +plt.plot(tout2, xout2[3,:], label=r'$k_{x, \gamma = 1.0}$') +plt.plot(tout3, xout3[3,:], label=r'$k_{x, \gamma = 5.0}$') +plt.hlines(kx_star, 0, Tend, label=r'$k_x^{\ast}$', color='black', linestyle='--') +plt.legend(loc=4, fontsize=14) +plt.title(r'control gain $k_x$ (feedback)') +plt.show() diff --git a/examples/mrac_siso_mit.py b/examples/mrac_siso_mit.py new file mode 100644 index 000000000..1f87dd5f6 --- /dev/null +++ b/examples/mrac_siso_mit.py @@ -0,0 +1,182 @@ +# mrac_siso_mit.py +# Johannes Kaisinger, 3 July 2023 +# +# Demonstrate a MRAC example for a SISO plant using MIT rule. +# Based on [1] Ex 5.2, Fig 5.5 & 5.6. +# Notation as in [2]. +# +# [1] K. J. Aström & B. Wittenmark "Adaptive Control" Second Edition, 2008. +# +# [2] Nhan T. Nguyen "Model-Reference Adaptive Control", 2018. + +import numpy as np +import scipy.signal as signal +import matplotlib.pyplot as plt + +import control as ct + +# Plant model as linear state-space system +A = -1. +B = 0.5 +C = 1 +D = 0 + +io_plant = ct.ss(A, B, C, D, + inputs=('u'), outputs=('x'), states=('x'), name='plant') + +# Reference model as linear state-space system +Am = -2 +Bm = 2 +Cm = 1 +Dm = 0 + +io_ref_model = ct.ss(Am, Bm, Cm, Dm, + inputs=('r'), outputs=('xm'), states=('xm'), name='ref_model') + +# Adaptive control law, u = kx*x + kr*r +kr_star = (Bm)/B +print(f"Optimal value for {kr_star = }") +kx_star = (Am-A)/B +print(f"Optimal value for {kx_star = }") + +def adaptive_controller_state(t, xc, uc, params): + """Internal state of adaptive controller, f(t,x,u;p)""" + + # Parameters + gam = params["gam"] + Am = params["Am"] + Bm = params["Bm"] + signB = params["signB"] + + # Controller inputs + r = uc[0] + xm = uc[1] + x = uc[2] + + # Controller states + x1 = xc[0] # + # x2 = xc[1] # kr + x3 = xc[2] # + # x4 = xc[3] # kx + + # Algebraic relationships + e = xm - x + + # Controller dynamics + d_x1 = Am*x1 + Am*r + d_x2 = - gam*x1*e*signB + d_x3 = Am*x3 + Am*x + d_x4 = - gam*x3*e*signB + + return [d_x1, d_x2, d_x3, d_x4] + +def adaptive_controller_output(t, xc, uc, params): + """Algebraic output from adaptive controller, g(t,x,u;p)""" + + # Controller inputs + r = uc[0] + # xm = uc[1] + x = uc[2] + + # Controller state + kr = xc[1] + kx = xc[3] + + # Control law + u = kx*x + kr*r + + return [u] + +params={"gam":1, "Am":Am, "Bm":Bm, "signB":np.sign(B)} + +io_controller = ct.nlsys( + adaptive_controller_state, + adaptive_controller_output, + inputs=('r', 'xm', 'x'), + outputs=('u'), + states=4, + params=params, + name='control', + dt=0 +) + +# Overall closed loop system +io_closed = ct.interconnect( + [io_plant, io_ref_model, io_controller], + connections=[ + ['plant.u', 'control.u'], + ['control.xm', 'ref_model.xm'], + ['control.x', 'plant.x'] + ], + inplist=['control.r', 'ref_model.r'], + outlist=['plant.x', 'control.u'], + dt=0 +) + +# Set simulation duration and time steps +Tend = 100 +dt = 0.1 + +# Define simulation time +t_vec = np.arange(0, Tend, dt) + +# Define control reference input +r_vec = np.zeros((2, len(t_vec))) +square = signal.square(2 * np.pi * 0.05 * t_vec) +r_vec[0, :] = square +r_vec[1, :] = r_vec[0, :] + +plt.figure(figsize=(16,8)) +plt.plot(t_vec, r_vec[0,:]) +plt.title(r'reference input $r$') +plt.show() + +# Set initial conditions, io_closed +X0 = np.zeros((6, 1)) +X0[0] = 0 # state of plant, (x) +X0[1] = 0 # state of ref_model, (xm) +X0[2] = 0 # state of controller, +X0[3] = 0 # state of controller, (kr) +X0[4] = 0 # state of controller, +X0[5] = 0 # state of controller, (kx) + +# Simulate the system with different gammas +tout1, yout1, xout1 = ct.input_output_response(io_closed, t_vec, r_vec, X0, + return_x=True, params={"gam":0.2}) +tout2, yout2, xout2 = ct.input_output_response(io_closed, t_vec, r_vec, X0, + return_x=True, params={"gam":1.0}) +tout3, yout3, xout3 = ct.input_output_response(io_closed, t_vec, r_vec, X0, + return_x=True, params={"gam":5.0}) + +plt.figure(figsize=(16,8)) +plt.subplot(2,1,1) +plt.plot(tout1, yout1[0,:], label=r'$x_{\gamma = 0.2}$') +plt.plot(tout2, yout2[0,:], label=r'$x_{\gamma = 1.0}$') +plt.plot(tout2, yout3[0,:], label=r'$x_{\gamma = 5.0}$') +plt.plot(tout1, xout1[1,:] ,label=r'$x_{m}$', color='black', linestyle='--') +plt.legend(fontsize=14) +plt.title(r'system response $x, (x_m)$') +plt.subplot(2,1,2) +plt.plot(tout1, yout1[1,:], label=r'$u_{\gamma = 0.2}$') +plt.plot(tout2, yout2[1,:], label=r'$u_{\gamma = 1.0}$') +plt.plot(tout3, yout3[1,:], label=r'$u_{\gamma = 5.0}$') +plt.legend(loc=4, fontsize=14) +plt.title(r'control $u$') +plt.show() + +plt.figure(figsize=(16,8)) +plt.subplot(2,1,1) +plt.plot(tout1, xout1[3,:], label=r'$k_{r, \gamma = 0.2}$') +plt.plot(tout2, xout2[3,:], label=r'$k_{r, \gamma = 1.0}$') +plt.plot(tout3, xout3[3,:], label=r'$k_{r, \gamma = 5.0}$') +plt.hlines(kr_star, 0, Tend, label=r'$k_r^{\ast}$', color='black', linestyle='--') +plt.legend(loc=4, fontsize=14) +plt.title(r'control gain $k_r$ (feedforward)') +plt.subplot(2,1,2) +plt.plot(tout1, xout1[5,:], label=r'$k_{x, \gamma = 0.2}$') +plt.plot(tout2, xout2[5,:], label=r'$k_{x, \gamma = 1.0}$') +plt.plot(tout3, xout3[5,:], label=r'$k_{x, \gamma = 5.0}$') +plt.hlines(kx_star, 0, Tend, label=r'$k_x^{\ast}$', color='black', linestyle='--') +plt.legend(loc=4, fontsize=14) +plt.title(r'control gain $k_x$ (feedback)') +plt.show() \ No newline at end of file