Figure 2.10:

Region of attraction (shaded region) for constrained MPC controller of Exercise \ref {exer:Xnterminal}.

Code for Figure 2.10

Text of the GNU GPL.

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
# [makes] Xnterminal_c.dat

# Exercise "Terminal constraint and region of attraction" part (c)
# Define model, cost function, and bounds.
import numpy as np
from misc import gnuplotsave
from findXn import findXn

def vertices2lines(p, xlims, ylims):
    """
    Returns 2-column array of line segments corresponding to edges of the
    polygon with vertices in columns of p (which must be in proper order).
    Each line segment has 101 interpolation points (to make gnuplot happy),
    and there is a row of NaNs inserted between each line.
    """
    dp = np.diff(p, axis=1)
    N = p.shape[1] - 1
    blocks = []
    xlims = np.linspace(xlims[0], xlims[1], 101).reshape(-1,1)
    ylims = np.linspace(ylims[0], ylims[1], 101).reshape(-1,1)

    for i in range(N):
        slope = dp[1,i]/dp[0,i]
        if np.isinf(slope) or np.isnan(slope):
            x = p[0,i] * np.ones((2,1))
            y = ylims
        else:
            x = xlims
            y = slope*(xlims - p[0,i]) + p[1,i]
        blocks.append(np.hstack((x, y)))
        blocks.append(np.array([[np.nan, np.nan]]))

    return np.vstack(blocks)

# Define model, cost function, and bounds
A = np.array([[2, 1], [0, 2]])
B = np.array([[1, 0], [0, 1]])
N = 3

alpha = 1e-5
Q = alpha * np.eye(2)
R = np.eye(2)

K = np.zeros(B.T.shape)

# Bounds
xlb = np.array([-np.inf, -np.inf])
xub = np.array([5, np.inf])
ulb = np.array([-1, -1])
uub = np.array([1, 1])

# Find region of stability
Xn, V, Z = findXn(A, B, K, N, xlb, xub, ulb, uub, 'termeq')

p = V[-1] # Vertices of Xn
xlims = [-3, 3]
ylims = [-3, 3]
Xnlines = vertices2lines(p, xlims, ylims)

# Save data files
p = p.T # Vertices of Xn
l = Xnlines # Lines for all sides

# Save to .dat file

gnuplotsave('Xnterminal_c.dat', {'p': p, 'l': l})