Figure 2.9:

Phase portrait for closed-loop evolution of cooler system with \dot {Q}_min=9. Line colors show value of discrete actuator u_2.

Code for Figure 2.9

Text of the GNU GPL.

getcooler.py

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
# sys = getcooler()
#
# Returns a struct of parameters for the cooler system.
import numpy as np
from scipy.linalg import expm, solve_discrete_lyapunov

def getcooler():
    """
    sys = getcooler()

    Returns a struct of parameters for the cooler system.
    """
    # Input validation
    if len(locals()) > 0:
        raise ValueError('Function takes no input arguments')

    # Define system
    Nx = 2
    Nu = 2
    alpha = 2
    beta = 1
    rho1 = 1
    rho2 = 1

    ns = np.array([0, 1, 2])
    Qmin = 9
    Qmax = 10

    nss = 1
    Qss = nss*Qmax

    uss = np.array([[Qss], [nss]])

    T0 = 40
    T1ss = 35
    T2ss = 25

    xss = np.array([[T1ss], [T2ss]])

    # Continuous-time system
    a = np.array([[-alpha - rho1, rho1],
                  [rho2, -rho2]])
    b = np.array([[0, 0],
                  [-beta, 0]])
    f = np.array([[alpha*T0],
                  [0]])

    # Discretize
    Delta = 0.25
    A = expm(Delta*a)
    B = np.linalg.solve(a, (A - np.eye(2)) @ b)
    F = np.linalg.solve(a, (A - np.eye(2)) @ f)

    # Choose penalty
    Q = np.eye(Nx)
    R = np.eye(Nu)

    # Pick terminal region
    Xf = {'rho': 1, 'P': solve_discrete_lyapunov(A.T, Q)}  # Level set for Xf

    # Save everything to a dict
    sys = {}
    local_vars = locals()
    for var in list(local_vars.keys()):
        if var != 'sys':
            sys[var] = local_vars[var]

    return sys

main.py

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
# Closed-loop optimization of cooler system.
# [makes] coolercl.mat
# [depends] functions/getcooler.py
## [rule] copymat
import sys
import os
sys.path.append('functions/')
import numpy as np
import scipy.io as sio
import matplotlib.pyplot as plt
import mpctools as mpc
from getcooler import getcooler

# get system
sys = getcooler()
Nx = sys['Nx']
Nu = sys['Nu']
xss = sys['xss']
uss = sys['uss']
F = sys['F'].flatten()  # avoid (2,) + (2,1) broadcast to (2,2) when mpctools
                        # wraps x, u as 1-D numpy arrays of CasADi scalars
f = mpc.getCasadiFunc(lambda x,u: sys['A'] @ x + sys['B'] @ u + F,
                      [Nx, Nu], ['x', 'u'], 'f')

# Choose cost function
Q = sys['Q']
R = 0.001 * sys['R']
P = sys['Xf']['P']
l = mpc.getCasadiFunc(lambda x, u, x_sp, u_sp: (x - x_sp).T @ Q @ (x - x_sp) + (u - u_sp).T @ R @ (u - u_sp),
                      [Nx, Nu, Nx, Nu], ['x', 'u', 'x_sp', 'u_sp'], 'l')

Vf = mpc.getCasadiFunc(lambda x, x_sp: (x - x_sp).T @ P @ (x - x_sp),
                       [Nx, Nx], ['x', 'x_sp'], 'Vf')

# Define input constraints and terminal constraint.
Qmin = sys['Qmin']
Qmax = sys['Qmax']
e = mpc.getCasadiFunc(lambda u: [u[0] - u[1] * Qmax,
                                 u[1] * Qmin - u[0]],
                      [Nu], ['u'], 'e')

udiscrete = np.array([False, True])

rho = sys['Xf']['rho'] # Level set for terminal region
ef = mpc.getCasadiFunc(lambda x, x_sp: (x - x_sp).T @ P @ (x - x_sp) - rho,
                       [Nx, Nx], ['x', 'x_sp'], 'ef')

# Build controller
Nt = 8
N = {}
N['x'] = Nx
N['u'] = Nu
N['t'] = Nt
N['e'] = 2
nmax = max(sys['ns'])

lb, ub = {}, {}
lb['u'] = np.zeros(Nu)
ub['u'] = np.array([nmax * Qmax, nmax])

init = {}
init['x_sp'] = xss
init['u_sp'] = uss

# mpc regulator
controller = mpc.nmpc(f=f,
                     l=l,
                     Vf=Vf,
                     N=N,
                     e=e,
                     ef=ef,
                     lb=lb, ub=ub,
                     par={'xsp': init['x_sp'], 'usp': init['u_sp']},
                     verbosity=0,
                     solver='bonmin',
                     udiscrete=udiscrete,
                     inferargs=True)
# Bonmin's DiveMIP heuristic asserts in some CasADi/Bonmin binary builds on
# small problems like this one; disable it so solve() doesn't abort.
controller.initialize(solveroptions={'heuristic_dive_MIP_fractional': 'no',
                                     'heuristic_dive_MIP_vectorLength': 'no',
                                     'heuristic_dive_fractional': 'no',
                                     'heuristic_dive_vectorLength': 'no'})

# Run closed-loop simulation for various initial conditions.
Nx0 = 16
theta = np.linspace(0, 2 * np.pi, Nx0 + 1)[1:]
x0 = xss + np.vstack([25 * np.cos(theta), 15 * np.sin(theta)])
Nsim = 16
x = np.full((Nx, Nsim + 1, Nx0), np.nan)
u = np.full((Nu, Nsim, Nx0), np.nan)

for j in range(Nx0):
    x[:, 0, j] = x0[:, j]
    for i in range(Nsim):
        controller.fixvar('x', 0, x[:, i, j])
        controller.solve()
        if controller.stats['status'] != 'SUCCESS':
            print(f'Solver failed at time {i} with status {controller.stats["status"]}')
            break
        u[:, i, j] = np.squeeze(controller.var['u', 0])
        x[:, i + 1, j] = np.squeeze(controller.var['x', 1])

lims = [sys['T1ss'] + 25 * np.array([-1, 1]), sys['T2ss'] + 15 * np.array([-1, 1])]
colors = ['r', 'c', 'b']
labels = ['n = 0', 'n = 1', 'n = 2']

if not os.getenv('OMIT_PLOTS') == 'true':
    plt.figure()
    plt.axis([lims[0][0], lims[0][1], lims[1][0], lims[1][1]])

    for j in range(Nx0):
        for i in range(Nsim):
            plt.plot(x[0, i:i+2, j], x[1, i:i+2, j], color=colors[int(u[1, i, j])])

    for k in range(3):
        plt.plot(np.nan, np.nan, '-' + colors[k], label=labels[k])

    plt.legend()
    plt.tight_layout()
    plt.show()