Figure 3.9:

Observed probability \varepsilon _test of constraint violation for i=10. Distribution is based on 500 trials for each value of \varepsilon . Dashed line shows the outcome predicted by formula \eqref {eq:tempoformula}, i.e., \varepsilon _test = \varepsilon .

Code for Figure 3.9

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
122
123
124
125
126
127
# Stochastic constraint tightening for a scalar system.
# [makes] sampling.mat
import numpy as np
import scipy.io
from scipy.linalg import solve_discrete_are
import matplotlib.pyplot as plt
from scipy.stats import norm, uniform
import os

# Stochastic constraint tightening for a scalar system.
# [makes] sampling.mat

A = np.array([[1]])
B = np.array([[1]])
Q = np.array([[0.5]])
R = np.array([[1]])

# Compute LQR gain
P = solve_discrete_are(A, B, Q, R)
K = -np.linalg.inv(R + B.T@P@B)@(B.T@P@A)
K = K[0,0] # Convert to scalar

AK = A[0,0] + B[0,0]*K

U = np.array([-1, 1])
X = np.array([-1, 1])
W = np.array([-0.5, 0.5])

N = 10  # Number of time points over which to consider disturbance

SKinf = W/(1 - AK)
Ubar = U - abs(K)*SKinf
Xbar_robust = X - SKinf

def evolution(A, B, Nt):
    """Computes the state evolution matrix for N steps."""
    A = np.atleast_2d(A)
    B = np.atleast_2d(B)
    Nx = B.shape[0]
    Nu = B.shape[1]
    E = np.zeros((Nt*Nx, Nt*Nu))
    e = np.hstack([B, np.zeros((Nx, (Nt-1)*Nu))])
    for k in range(Nt):
        i = slice(k*Nx, (k+1)*Nx)
        E[i,:] = e
        if k < Nt-1:
            e = A@e
            j = slice((k+1)*Nu, (k+2)*Nu)
            e[:,j] = B
    return E

E = evolution(AK, 1, N)  # Converts from vector of w(k) to vector of e(k)
np.random.seed(0)

def sampleE(n):
    return E@(W[0] + (W[1] - W[0])*np.random.rand(N, n))

# Reduce number of trials and samples for performance
Ntrials = 500
beta = 0.01
epsilons = np.logspace(-3, 0, 16)
Ms = np.ceil(2./epsilons*(2 - np.log(beta))).astype(int)

esim = sampleE(1000000).reshape(-1, 1000000)

Xbar = np.full((2, len(Ms), Ntrials), np.nan)
eviolation = np.full((len(Ms), Ntrials), np.nan)

for i in range(len(Ms)):
    print(f'{i:2d}. epsilon = {epsilons[i]:10g} (M = {Ms[i]})')
    for j in range(Ntrials):
        e = sampleE(Ms[i]).reshape(-1, Ms[i])
        emin = min(np.min(e[0,:]), 0)
        emax = max(np.max(e[0,:]), 0)
        eviolation[i,j] = np.mean(esim[0,:] < emin) + np.mean(esim[0,:] > emax)
        Xbar[0,i,j] = min(X[0] - emin, 0)
        Xbar[1,i,j] = max(X[1] - emax, 0)

def violationquantiles(x, v, logscale=True):
    """Plots quantiles of violation probability."""
    quantiles = np.array([0.05, 0.25, 0.5, 0.75, 0.95])
    colors = plt.cm.jet(np.linspace(0,1,len(quantiles)))
    q = np.quantile(np.squeeze(v), quantiles, axis=1).T

    if not os.getenv('OMIT_PLOTS') == 'true':
        if logscale:
            q = np.maximum(q, 1e-16)
            plotfunc = plt.loglog
        else:
            plotfunc = plt.semilogx

        plt.figure()
        for j in range(len(quantiles)):
            plotfunc(x, q[:,j], '-o', color=colors[j],
                    label=f'percentile {100*quantiles[j]:.0f}')
        plt.legend()
    return q

if not os.getenv('OMIT_PLOTS') == 'true':
    # Plot tightened sets
    labels = ['Xbar_{lb}', 'Xbar_{ub}']
    locations = ['lower left', 'upper left']
    for i in range(2):
        plt.figure()
        q = violationquantiles(epsilons, Xbar[i,:,:], False)
        plt.plot(epsilons[[0,-1]], Xbar_robust[i]*np.ones(2), '-k', label='Robust')
        plt.xlabel('ε')
        plt.ylabel(labels[i])
        plt.legend(loc=locations[i])
        plt.title('Tightened Constraints')
        plt.axis([min(epsilons), max(epsilons), i-2.1, i-0.9])

    # Plot violation probabilities
    plt.figure()
    q = violationquantiles(epsilons, eviolation, True)
    plt.plot(epsilons, epsilons, '--k')
    plt.legend(loc='lower right')
    plt.ylabel('ε_{test}')
    plt.xlabel('ε')
    plt.axis([min(epsilons), max(epsilons), 1e-6, 1])

    plt.show()

# Save data
data = {'Ms': Ms, 'Xbar': Xbar, 'Xbar_robust': Xbar_robust,
        'eviolation': eviolation, 'beta': beta, 'epsilons': epsilons}
scipy.io.savemat('sampling.mat', data)