Figure 2.2:

Feasible region \mathcal {U}_2, elliptical cost contours and ellipse center a(x), and constrained minimizers for different values of x.

Code for Figure 2.2

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
67
68
69
70
71
72
73
74
75
76
77
78
# Exercise for ch2; Continuation of Example 2.5: Linear quadratic MPC
# jbr, March 1, 2009

import numpy as np
from scipy import linalg
import matplotlib.pyplot as plt
from misc import ellipse
import os

H = np.array([[3, 1], [1, 2]])
hr1 = H[0,:]
hr2 = H[1,:]

Hinv = (1/5)*np.array([[2, -1], [-1, 3]])
n1 = 1/np.sqrt(5)*np.array([2, -1])

lam, V = linalg.eig(H)
v1 = V[:,0]
v2 = V[:,1]

box = np.array([[-1, -1], [1, -1], [1, 1], [-1, 1], [-1, -1]])

umin = -4.75
av = -np.array([3, 1])
line = np.vstack((0*av, -umin/3*av))

xc1 = 5/3

xc2 = 3
ac2 = xc2/5*av
bc2 = 5/2*(3*xc2/5-1)**2
xec2, yec2, *_ = ellipse(H, bc2, 100, ac2)

x3 = 4.5
a3 = x3/5*av
b3 = (np.array([-1,-1])-a3).T @ H @ (np.array([-1,-1])-a3)
xe3, ye3, *_ = ellipse(H, b3, 100, a3)

x4 = 2.25
a4 = x4/5*av
b4 = 5/2*(3*x4/5-1)**2
xe4, ye4, *_ = ellipse(H, b4, 100, a4)
u0 = np.array([-1, -(x4-1)/2])

if not os.getenv('OMIT_PLOTS') == 'true':
    plt.figure()
    plt.plot(box[:,0], box[:,1], '-',
             0, 0, 'o',
             -1, -1/3, 'o',
             line[:,0], line[:,1], '-',
             ac2[0], ac2[1], 'x',
             xec2, yec2, '-',
             a3[0], a3[1], 'x',
             xe3, ye3, '-',
             -1, -1, 'o',
             a4[0], a4[1], 'x',
             xe4, ye4, '-',
             u0[0], u0[1], 'o')

    plt.axis([umin, 1.5, -3.5, 1.5])
    plt.show(block=False)

table1 = np.column_stack((xec2, yec2, xe3, ye3, xe4, ye4))
table2 = np.hstack(([0, 0], [-1, -1/3], ac2, a3, a4, [-1, -1], u0))

# Save data
with open('constlq.dat', 'w') as f:
    f.write('# table1\n')
    np.savetxt(f, table1, delimiter=' ', fmt='%.6f')
    f.write('\n\n')
    f.write('# table2\n')
    np.savetxt(f, [table2], delimiter=' ', fmt='%.6f')
    f.write('\n\n')
    f.write('# box\n')
    np.savetxt(f, box, delimiter=' ', fmt='%.6f')
    f.write('\n\n')
    f.write('# line\n')
    np.savetxt(f, line, delimiter=' ', fmt='%.6f')