Figure 6.5:

Ten iterations of cooperative steady-state calculation; reversed pairing.

Code for Figure 6.5

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
# A 2-variable steady-state problem.
import os
import matplotlib.pyplot as plt
import numpy as np
from scipy import linalg

from misc import gnuplotsave

# A 2-variable steady-state problem
lams = [[0.5, 0.5], [0.2, 0.8], [0.8, 0.2], [0.5, 0.5]]
iters = [10, 25, 25, 25]
names = ['equal', 'first', 'second', 'swapped']
swaps = [False, False, False, True]

def noncoopvscoop(lam, Niter, name, swap=False):
    """Simulates non-cooperative vs. cooperative optimization."""

    # Unstable for about m1 < 0.58, lam1 = lam2 = 0.5
    m1 = 0.5
    m2 = -1/m1
    g = np.array([[-m1, 1], [-m2, 1]])
    if swap:
        g = np.array([[0, 1], [1, 0]]) @ g

    ysp = np.array([[1], [1]])
    usp = linalg.solve(g, ysp)
    ude = np.array([[ysp[0,0]/g[0,0]], [ysp[1,0]/g[1,1]]])

    # Noncooperative problem
    Lnc = np.array([[1-lam[0], -lam[0]*g[0,1]/g[0,0]],
                    [-lam[1]*g[1,0]/g[1,1], 1-lam[1]]])
    Knc = np.diag([lam[0]/g[0,0], lam[1]/g[1,1]])

    ssnc = np.abs(linalg.eigvals(Lnc))

    unc = np.zeros((2, Niter + 1))
    unc[:,0] = ude.flatten()

    for i in range(Niter):
        unc[:,i+1] = (Knc @ ysp + Lnc @ unc[:,i].reshape(-1,1)).flatten()

    # Cooperative problem
    Q = np.diag(lam)
    H = np.block([[Q, np.zeros((2,1))],
                  [np.zeros((1,2)), 0]])
    G1 = g[:,0:1]
    G2 = g[:,1:2]
    D = np.hstack([np.eye(2), -G1])
    M = np.block([[H, -D.T], [-D, np.zeros((2,2))]])
    Minv = linalg.inv(M)
    K1 = Minv[2,0:2] @ Q
    L1 = -Minv[2,3:5] @ G2
    D = np.hstack([np.eye(2), -G2])
    M = np.block([[H, -D.T], [-D, np.zeros((2,2))]])
    Minv = linalg.inv(M)
    K2 = Minv[2,0:2] @ Q
    L2 = -Minv[2,3:5] @ G1
    Kco = np.vstack([lam[0]*K1, lam[1]*K2])
    Lco = np.array([[1-lam[0], lam[0]*L1[0]],
                    [lam[1]*L2[0], 1-lam[1]]])
    ssco = np.abs(linalg.eigvals(Lco))

    uco = np.zeros((2, Niter+1))
    uco[:,0] = ude.flatten()

    for i in range(Niter):
        uco[:,i+1] = (Kco @ ysp + Lco @ uco[:,i].reshape(-1,1)).flatten()

    # Make lines for plotting
    ld = np.min(unc, axis=1)
    ru = np.max(unc, axis=1)
    lr = max(-ld[0], ru[0])
    ud = max(-ld[1], ru[1])
    sslines = np.array([[-lr, lr, np.nan, -ud/m2, ud/m2],
                       [-m1*lr, m1*lr, np.nan, -ud, ud]])
    sslines = sslines + np.tile(usp, (1,5))

    if not os.getenv('OMIT_PLOTS') == 'true':
        plt.figure()
        plt.plot(sslines[0,:], sslines[1,:], '-k',
                 unc[0,:], unc[1,:], '-sr',
                 uco[0,:], uco[1,:], '-og')
        plt.title(f'{name} (lambda = [{lam[0]}, {lam[1]}])')
        plt.legend(['Steady State', 'Non-cooperative', 'Cooperative'])
        plt.xlabel('u_1')
        plt.ylabel('u_2')

    # Package data
    data = {'ss': sslines.T, 'noncoop': unc.T, 'coop': uco.T}
    return data

# Simulate with different values of lambda
data = {}
for i in range(len(lams)):
    data[names[i]] = noncoopvscoop(lams[i], iters[i], names[i], swaps[i])

# Save data
gnuplotsave('stst.dat', data)