import numpy as np
import cvxpy as cvx
import pandas as pd
from matplotlib import pyplot as plt
n = 32
m = 8
P = np.array([ [ 1/m for i in range(n) ] for j in range(m) ])
b = np.array([1e3 for i in range(m)])
c_hat = np.array([0.8 + 0.2*np.sqrt((j-1)/(n-1)) for j in np.arange(1,n+1)])
alphas = np.array([0, 0.001, 0.01, 0.1, 1, 10])
def compute(n, m, P, b, c_hat, alpha):
"""Solve the optimization problem.
"""
sigma = np.array([alpha*(1.2 - c_hat[j]) for j in range(n)])
X = cvx.Variable((n))
z = cvx.Variable((n))
obj = cvx.Minimize(
cvx.sum(c_hat @ X + z)
)
constraints = [
P @ X >= b,
X >= 0,
z >= cvx.norm(cvx.diag(sigma @ X), 2)
]
prob = cvx.Problem(obj, constraints)
min_cost = prob.solve()
return X, min_cost
amounts, min_cost = [], []
for alpha in alphas:
x_arr, minC = compute(n, m, P, b, c_hat, alpha)
amounts.append(x_arr.value)
min_cost.append(minC)
print(min_cost)
print(alphas)
[204806.67568375397, 204906.67136196772, 205823.24017668262, 215035.18805363585, 307199.98716804455, 768010.3296590501]
[0.e+00 1.e-03 1.e-02 1.e-01 1.e+00 1.e+01]
plt.plot(alphas, min_cost)
plt.xlabel("Volatility (alpha)")
plt.ylabel("Optimal Supply Costs")
plt.title('Trade-Off Curve');