User Guide for Python+RSOME

整理自 RSOME for Python：所有例子都由熊老师的网站提供，本笔记只是整理一下方便阅读（暂时没有PDF格式的User Guide）。

1. Deterministic linear and second-order cone programs

1.1 Mean-variance portfolio optimization

import rsome as rso
import numpy as np
from rsome import ro
from rsome import grb_solver as grb

n = 150                                     # Number of stocks
i = np.arange(1, n+1)                       # Indices of stocks
p = 1.15 + i*0.05/150                       # Mean returns
sigma = 0.05/450 * (2*i*n*(n+1))**0.5       # Standard deviations of returns
phi = 5                                     # Constant phi

model = ro.Model('mv-portfolio')

x = model.dvar(n)                           # Fractions of investment

model.max(p@x - phi*rso.sumsqr(sigma*x))    # Mean-variance objective
model.st(x.sum() == 1)                      # Summation of x is one
model.st(x >= 0)                            # x is non-negative

model.solve(grb)

import matplotlib.pyplot as plt

obj_val = model.get()               # The optimal objective value
x_sol = x.get()                     # The optimal investment decision

plt.plot(range(1, n+1), x_sol,
         linewidth=2, color='b')
plt.xlabel('Stocks')
plt.ylabel('Fraction of investment')
plt.show()
print('Objective value: {0:0.4f}'.format(obj_val))

1.2 Integer programming for Sudoku

import rsome as rso
import numpy as np
from rsome import ro
from rsome import grb_solver as grb

# A Sudoku puzzle
# Zeros represent unknown numbers
puzzle = np.array([[5, 3, 0, 0, 7, 0, 0, 0, 2],
                   [6, 0, 0, 1, 9, 5, 0, 0, 0],
                   [0, 9, 8, 0, 0, 0, 0, 6, 0],
                   [8, 0, 0, 0, 6, 0, 0, 0, 3],
                   [4, 0, 0, 8, 0, 3, 0, 0, 1],
                   [7, 0, 0, 0, 2, 0, 0, 0, 6],
                   [0, 6, 0, 0, 0, 0, 2, 8, 0],
                   [0, 0, 0, 4, 1, 9, 0, 0, 5],
                   [0, 0, 0, 0, 8, 0, 0, 7, 9]])

# Create model and binary decision variables
model = ro.Model()
x = model.dvar((9, 9, 9), vtype='B')

# Objective is set to be zero
model.min(0 * x.sum())

# Constraints 1 to 3
model.st(x.sum(axis=0) == 1,
         x.sum(axis=1) == 1,
         x.sum(axis=2) == 1)

# Constraints 4
i, j = np.where(puzzle)
model.st(x[i, j, puzzle[i, j]-1] == 1)

# Constraints 5
for i in range(0, 9, 3):
    for j in range(0, 9, 3):
        model.st(x[i: i+3, j: j+3, :].sum(axis=(0, 1)) == 1)

# Solve the integer programming problem
model.solve(grb)

2. Robust optimizaiton models

2.1 Robust portfolio optimization

from rsome import ro
from rsome import grb_solver as grb
import rsome as rso
import numpy as np


n = 150                                 # Number of stocks
i = np.arange(1, n+1)                   # Indices of stocks
p = 1.15 + i*0.05/150                   # Mean returns
delta = 0.05/450 * (2*i*n*(n+1))**0.5   # Deviations of returns
Gamma = 5                               # Budget of uncertainty

model = ro.Model()              
x = model.dvar(n)                       # Fractions of investment
z = model.rvar(n)                       # Random variables

model.maxmin((p + delta*z) @ x,         # The max-min objective
             rso.norm(z, np.infty) <=1, # Uncertainty set constraints
             rso.norm(z, 1) <= Gamma)   # Uncertainty set constraints
model.st(sum(x) == 1)                   # Summation of x is one
model.st(x >= 0)                        # x is non-negative

model.solve(grb)                        # Solve the model by Gurobi

2.2 The robust and robustness knapsack problems

from rsome import ro
from rsome import grb_solver as grb
import rsome as rso
import numpy as np
import numpy.random as rd
import matplotlib.pyplot as plt

N = 50
b = 2000

c = 2*rd.randint(low=5, high=10, size=N)    # Profit coefficients
w = 2*rd.randint(low=10, high=41, size=N)   # Nominal weights
delta = 0.2*w                               # Maximum deviations

def robust(r):
    """
    The function robust implements the robust optimiztion model,
    given the budget of uncertainty r
    """

    model = ro.Model('robust')
    x = model.dvar(N, vtype='B')    
    z = model.rvar(N)              

    z_set = (abs(z) <= 1, rso.norm(z, 1) <= r)
    model.max(c @ x)
    model.st(((w + z*delta) @ x <= b).forall(z_set))

    model.solve(grb, display=False) # Disable solution message

    return model.get(), x.get()     # Optimal objective and solution

def robustness(tau):
    """
    The function robustness implements the robustness optimization
    model, given the profit target tau.
    """

    model = ro.Model('robustness')

    x = model.dvar(N, vtype='B')    
    k = model.dvar()              
    z = model.rvar(N)           
    u = model.rvar(N)           

    z_set = (abs(z) <= u, u <= 1)
    model.min(k)
    model.st(c @ x >= tau)
    model.st(((w + z*delta) @ x - b <= k*u.sum()).forall(z_set))
    model.st(k >= 0)

    model.solve(grb, display=False) # Disable solution message

    return model.get(), x.get()     # Optimal objective and solution

def sim(x_sol, zs):
    """
    The function sim is for calculating the probability of violation
    via simulations.
        x_sol: solution of the Knapsack problem
        zs: random sample of the random variable z
    """

    ws = w + zs*delta   # Random samples of uncertain weights

    return (ws @ x_sol > b).mean()

step = 0.1
rs = np.arange(1, 5+step, step)         # All budgets of uncertainty
num_samp = 20000
zs = 1-2*rd.rand(num_samp, N)           # Random samples for z

"""Robust optimization"""
outputs_rb = [robust(r) for r in rs]
targets = [output[0]
           for output in outputs_rb]    # RO Objective as targets
pv_rb = [sim(output[1], zs)
         for output in outputs_rb]      # Prob. of violations

"""Robustness optimization"""
outputs_rbn = [robustness(target)
               for target in targets]   
pv_rbn = [sim(output[1], zs)
          for output in outputs_rbn]    # Prob. of violations

step = 0.1
rs = np.arange(1, 5+step, step)         # All budgets of uncertainty
num_samp = 20000
zs = 1-2*rd.rand(num_samp, N)           # Random samples for z

"""Robust optimization"""
outputs_rb = [robust(r) for r in rs]
targets = [output[0]
           for output in outputs_rb]    # RO Objective as targets
pv_rb = [sim(output[1], zs)
         for output in outputs_rb]      # Prob. of violations

"""Robustness optimization"""
outputs_rbn = [robustness(target)
               for target in targets]   
pv_rbn = [sim(output[1], zs)
          for output in outputs_rbn]    # Prob. of violations

plt.plot(rs, pv_rb, marker='o', markersize=5, c='b',
         label='Robust Optimization')
plt.plot(rs, pv_rbn, c='r',
         label='Robustness Optimization')

plt.legend()
plt.xlabel('Parameter r in robust optimization')
plt.ylabel('Prob. violation')
plt.show()

plt.scatter(targets, pv_rb, c='b', alpha=0.3,
            label='Robust Optimization')
plt.scatter(targets, pv_rbn, c='r', alpha=0.3,
            label='Robustness Optimization')

plt.legend()
plt.xlabel(r'Target return $\tau$')
plt.ylabel('Prob. violation')
plt.show()

2.3 Adaptive robust optimization for a lot-sizing problem

from rsome import ro
import rsome.grb_solver as grb
import numpy as np
import numpy.random as rd
import matplotlib.pyplot as plt

# Parameters of the lot-sizing problem
N = 30
c = 20
dmax = 20
Gamma = 20*np.sqrt(N)
xy = 10*rd.rand(2, N)
t = ((xy[[0]] - xy[[0]].T) ** 2
     + (xy[[1]] - xy[[1]].T) ** 2) ** 0.5

model = ro.Model()      # Define a model
x = model.dvar(N)       # Define location decisions
y = model.ldr((N, N))   # Define decision rule as the shifted quantities
d = model.rvar(N)       # Define random variables as the demand

y.adapt(d)              # The decision rule y affinely depends on d

# Define the uncertainty set
uset = (d >= 0, d <= dmax, sum(d) <= Gamma)

# Define model objective and constraints
model.minmax((c*x).sum() + (t*y).sum(), uset)
model.st(d <= y.sum(axis=0) - y.sum(axis=1) + x)
model.st(y >= 0)
model.st(x >= 0)
model.st(x <= 20)

model.solve(grb)

plt.figure(figsize=(5, 5))
plt.scatter(xy[0], xy[1], c='w', edgecolors='k')
plt.scatter(xy[0], xy[1], s=40*x.get(), c='k', alpha=0.6)
plt.axis('equal')
plt.xlim([-1, 11])
plt.ylim([-1, 11])
plt.grid()
plt.show()

3. Distributionally robust optimizaiton models

3.1 Distributionally robust optimization for medical appointment scheduling

from rsome import dro
from rsome import square
from rsome import E
from rsome import grb_solver as grb
import numpy as np
import numpy.random as rd


# Model and ambiguity set parameters
N = 8
gamma = 2
alpha = 0.25
mu = 30 + 30*rd.rand(N)
sigma = mu * 0.3*rd.rand(1, N)
T = mu.sum() + 0.5*((sigma**2).sum())**0.5

mul = np.diag(np.ones(N))*(1-alpha) + np.ones((N, N))*alpha
Sigma = (sigma.T @ sigma) * mul

# Modeling with RSOMEb
model = dro.Model()                                  # Define a DRO model
z = model.rvar(N)                                    # Random variable z
u = model.rvar(N+1)                                  # Auxiliary variable u

fset = model.wks(z >= 0, square(z - mu) <= u[:-1],
                 square((z-mu).sum()) <= u[-1],      # Support of random variables
                 E(z) == mu, E(u[:-1]) <= sigma**2,
                 E(u[-1]) <= Sigma.sum())            # Support of expectations

x = model.dvar(N)                                    # The first-stage decision
y = model.dvar(N+1)                                  # The decision rule
y.adapt(z)                                           # Define affine adaptation
y.adapt(u)                                           # Define affine adaptation

model.minsup(E(y[:-1].sum() + gamma*y[-1]), fset)    # Worst-case expected cost
model.st(y[1:] - y[:-1] + x >= z)                    # Constraints
model.st(y >= 0)                                     # Constraints
model.st(x >= 0)                                     # Constraints
model.st(x.sum() <= T)                               # Constraints

model.solve(grb)                                     # Solve the model by Gurobi

3.2 A multi-item newsvendor problem considering the Wasserstein ambiguity set

from rsome import dro
from rsome import norm
from rsome import E
from rsome import grb_solver as grb
import numpy as np
import numpy.random as rd

# Model and ambiguity set parameters
I = 2
S = 50
c = np.ones(I)
d = 50 * I
p = 1 + 4*rd.rand(I)
zbar = 100 * rd.rand(I)
zhat = zbar * rd.rand(S, I)
theta = 0.01 * zbar.min()

# Modeling with RSOME
model = dro.Model(S)                        # Create a DRO model with S scenarios
z = model.rvar(I)                           # Random demand z
u = model.rvar()                            # Auxiliary random variable

fset = model.ambiguity()                    #
for s in range(S):
    fset[s].suppset(0 <= z, z <= zbar,
                    norm(z - zhat[s]) <= u) # Define the support for each scenario
fset.exptset(E(u) <= theta)                 # The Wasserstein metric constraint
pr = model.p                                # An array of scenario probabilities
fset.probset(pr == 1/S)                     # Support of scenario probabilities

x = model.dvar(I)                           # Define first-stage decisions
y = model.dvar(I)                           # Define decision rule variables
y.adapt(z)                                  # y affinely adapts to z
y.adapt(u)                                  # y affinely adapts to u
for s in range(S):
    y.adapt(s)                              # y adapts to each scenario s

model.minsup(-p@x + E(p@y), fset)           # Worst-case expectation over fset
model.st(y >= 0)                            # Constraints
model.st(y >= x - z)                        # Constraints
model.st(x >= 0)                            # Constraints
model.st(c@x == d)                          # Constraints

model.solve(grb)                            # Solve the model by Gurobi