Figure 8.7:

Relaxed and binary feasible solution for Example~\ref {ex:nc-binary}.

Code for Figure 8.7

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
 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
# [makes] CIA.mat
## [rule] copymat
# Solve a small MINLP problem with a CIA (Combinatorial Integral Approximation)
# heuristic. System: x^+ = RK4(x, u, Delta) with f(x,u) = x^3 - u.
import numpy as np
from scipy.optimize import minimize, milp, LinearConstraint, Bounds
from scipy.io import savemat

# Problem parameters.
xref        = 0.7
Delta       = 0.05
x0          = 0.9
Ncontroller = 30
time        = Delta * np.arange(Ncontroller)


def rk4_step(x, u, h):
    """One step of RK4 for dx/dt = x^3 - u."""
    f = lambda xv: xv**3 - u
    k1 = f(x)
    k2 = f(x + h/2*k1)
    k3 = f(x + h/2*k2)
    k4 = f(x + h*k3)
    return x + h/6*(k1 + 2*k2 + 2*k3 + k4)


def simulate(u, x_init):
    """Simulate N steps, return state trajectory (N+1,)."""
    x = np.empty(len(u) + 1)
    x[0] = x_init
    for k in range(len(u)):
        x[k+1] = rk4_step(x[k], u[k], Delta)
    return x


# ── NLP: relax u to [0,1] and solve as a multiple-shooting problem ───────────
# Variables: z = [u_0,...,u_{N-1}, x_1,...,x_N] (length 2*N)
# Constraints: x_{k+1} - rk4_step(x_k, u_k, Delta) = 0 (N equalities)
Nu = Ncontroller
Nz = 2 * Nu

def unpack(z):
    return z[:Nu], z[Nu:]          # u (Nu,), x (Nu,) representing x_1..x_N

def nlp_obj_ms(z):
    u, x = unpack(z)
    x_full = np.concatenate(([x0], x))   # x_0..x_N
    return float(np.sum((x_full - xref)**2))

def nlp_defect(z):
    u, x = unpack(z)
    x_prev = np.concatenate(([x0], x[:-1]))  # x_k for k=0..N-1
    return x - np.array([rk4_step(xp, uk, Delta) for xp, uk in zip(x_prev, u)])

u_init_ms = np.ones(Nu)
x_init_ms = simulate(u_init_ms, x0)[1:]    # trajectory from constant u=1
z_init    = np.concatenate([u_init_ms, x_init_ms])

bounds_ms = [(0.0, 1.0)] * Nu + [(None, None)] * Nu
nlp_res   = minimize(nlp_obj_ms, z_init, method='SLSQP', bounds=bounds_ms,
                     constraints={'type': 'eq', 'fun': nlp_defect},
                     options={'ftol': 1e-12, 'maxiter': 500})
unlp = nlp_res.x[:Nu]
xnlp = simulate(unlp, x0)
Vnlp = nlp_res.fun


# ── MILP: CIA heuristic ───────────────────────────────────────────────────────
Nc = Ncontroller

A_lower = np.tril(np.ones((Nc, Nc)))  # cumulative-sum matrix
C_vec   = np.ones(Nc)

# B encodes switching constraints (from MATLAB code).
B = np.zeros((Nc, Nc))
for i in range(Nc - 2):
    B[i:i+3, i] = [-1, 1, -1]
B[Nc-2:Nc, Nc-2] = [-1, 1]
B[Nc-1,    Nc-1] = -1

# Variables: x = [eta, u_0, ..., u_{Nc-1}]
# Minimise eta.
c_obj = np.zeros(Nc + 1)
c_obj[0] = 1.0

# Constraints (all of the form lb <= A*x <= ub):
#   [-1, A_lower] * [eta;u] <=  A_lower @ unlp    (A*(u-unlp) <= eta)
#   [-1,-A_lower] * [eta;u] <= -A_lower @ unlp    (A*(unlp-u) <= eta)
#   [0,       B ] * [eta;u] <=  0                 (B*u <= 0)
A_con = np.zeros((3*Nc, Nc + 1))
A_con[:Nc,      0]  = -1.0;  A_con[:Nc,     1:] =  A_lower
A_con[Nc:2*Nc,  0]  = -1.0;  A_con[Nc:2*Nc, 1:] = -A_lower
A_con[2*Nc:,    1:]  = B

b_ub = np.concatenate([A_lower @ unlp, -(A_lower @ unlp), np.zeros(Nc)])

milp_bounds = Bounds(
    lb=np.concatenate([[-np.inf], np.zeros(Nc)]),
    ub=np.concatenate([[ np.inf], np.ones(Nc)]),
)
integrality = np.zeros(Nc + 1)
integrality[1:] = 1   # u is binary

milp_res     = milp(c_obj,
                    constraints=LinearConstraint(A_con, -np.inf, b_ub),
                    integrality=integrality,
                    bounds=milp_bounds)
maxcianorm = milp_res.x[0]
umilp      = milp_res.x[1:]

# Simulate with binary control.
xopt = simulate(umilp, x0)

savemat('CIA.mat', {
    'xref': xref,
    'time': time,
    'xnlp': xnlp,
    'unlp': unlp,
    'xopt': xopt,
    'umilp': umilp,
})