Production Optimization

Consider a simple production optimization problem where the optimal production amount \(x\) is to be determined subject to uncertainty in the demand \(\mathbf{d}\) which is represented by a random variable.

The optimization problem is given by:

\[\begin{split}\begin{align*} \min_x \quad & c \cdot x \\ \text{subject to} \quad & \mathbb{P}(x \geq \mathbf{d}) \geq 1 - \epsilon \end{align*}\end{split}\]

Where \(c\) is the unit cost of the product, \(\mathbb{P}\) measures the probability of an event, and \(\epsilon \in (0, 0.5]\) is the risk bound.

Assume that \(\mathbf{d}\) follows a Gaussian distribution with mean \(\overline{d}\) and variance \(\sigma^2\): \(\mathbf{d} \sim \mathcal{N}(\overline{d}, \sigma^2)\).

The chance constraint can be reformulated into a deterministic constraint as follows:

\[\begin{split}\begin{align*} \mathbb{P}(x \geq \mathbf{d}) \geq 1 - \epsilon & \iff \mathbb{P}(-x + \mathbf{d} \leq 0) \geq 1 - \epsilon \\ & \iff -x + \overline{d} + \sigma \Phi^{-1}(1-\epsilon) \leq 0 \\ & \iff x \geq \overline{d} + \sigma \Phi^{-1}(1-\epsilon) \end{align*}\end{split}\]

where \(\Phi^{-1}\) is the inverse CDF of the standard normal Gaussian distribution.

Example

In the following code, we solve the production optimization problem with CVXPY and cvxRiskOpt and with CVXPY only. The main difference between using cvxRiskOpt and not doing so is in the inclusion of the chance constraint.

With cvxRiskOpt we only need to rearrange the chance constraint:

\[\begin{split}\begin{align*} & \mathbb{P}(-x + \mathbf{d} \leq 0) \geq 1 - \epsilon \iff \mathbb{P}(a \mathbf{\xi_1} + b + \mathbf{\xi_2} \leq 0) \geq 1 - \epsilon\\ \rightarrow & a = 1, \ \mathbf{\xi_1} = \mathbf{d}, \ b=-x, \ \mathbf{\xi_2} = 0 \end{align*}\end{split}\]

However, using CVXPY only, we need to reformulate it into the deterministic constraint as show earlier.

import cvxpy as cp
import numpy as np
import time

solver = cp.OSQP

c_val = 10  # cost
x = cp.Variable(name='x')  # decision variable
c = cp.Parameter(name='c')  # create a parameter for the cost
d_mean = 700  # demain mean
d_var = 30  # demand variance
eps = 0.1  # risk bound

# cvxpy problems
objective = cp.Minimize(c * x)

# cvxpy + cvxRiskOpt
from cvxRiskOpt.cclp_risk_opt import cclp_gauss
cc_contr = cclp_gauss(eps, a=1, b=-x,
                      xi1_hat=d_mean,
                      gam11=d_var)
constraints_with_cro = [x >= 0, cc_contr]
prob_with_cro = cp.Problem(objective, constraints_with_cro)
c.value = c_val
prob_with_cro.solve(solver=solver)  # cvxpy's first run is usually slower than all other solves. Solve once before timing it
t0 = time.time()
prob_with_cro.solve(solver=solver)
t1 = time.time()
print("Production amount (CVXPY + cvxRiskOpt): ", x.value)
print("Solve Time: %.3f ms" % (1000 * (t1 - t0)))

# cvxpy only
from scipy.stats import norm
d_std_div = np.sqrt(d_var)
constraints = [x >= 0, x >= d_mean + d_std_div * norm.ppf(1-eps)]
prob = cp.Problem(objective, constraints)
c.value = c_val
prob.solve(solver=solver)  # cvxpy's first run is usually slower than all other solves. Solve once before timing it
t0 = time.time()
prob.solve(solver=solver)
t1 = time.time()
print("Production amount (CVXPY only): ", x.value)
print("Solve Time: %.3f ms" % (1000 * (t1 - t0)))

Production amount (CVXPY + cvxRiskOpt):  707.0193470105485
Solve Time: 0.572 ms
Production amount (CVXPY only):  707.0193470105482
Solve Time: 0.563 ms

Another benefit or using cvxRiskOpt is that if the Guassian assumption about the noise were to change, updating the chance constraint can easily be done by changing the cclp_gauss call to call another cclp_risk_opt function. Below is an example using a DR-VaR risk metric where the probability of meeting the demand must be realized under the work case distribution in a moment-based ambiguity set using the mean and covariance of the uncertainty

from cvxRiskOpt.cclp_risk_opt import cclp_dro_mean_cov
dr_cc_contr = cclp_dro_mean_cov(eps, a=1, b=-x,
                                xi1_hat=d_mean,
                                gam11=d_var)
dr_constraints_with_cro = [x >= 0, dr_cc_contr]
dr_prob_with_cro = cp.Problem(objective, dr_constraints_with_cro)
c.value = c_val
dr_prob_with_cro.solve(solver=solver)  # cvxpy's first run is usually slower than all other solves. Solve once before timing it
t0 = time.time()
dr_prob_with_cro.solve(solver=solver)
t1 = time.time()
print("DR Production amount (CVXPY + cvxRiskOpt): ", x.value)
print("Solve Time: %.3f ms" % (1000 * (t1 - t0)))

DR Production amount (CVXPY + cvxRiskOpt):  716.4316767251547
Solve Time: 0.603 ms

Generating C code

We can also generate C code for the CVXPY Problem instance above using CVXPYgen. This can be done as shown below.

from cvxpygen import cpg
cpg.generate_code(prob_with_cro, code_dir='prod_opt', solver=solver)
from prod_opt.cpg_solver import cpg_solve
prob_with_cro.register_solve('cpg', cpg_solve)
c.value = c_val
t0 = time.time()
prob_with_cro.solve(method='cpg', updated_params=['c'])
t1 = time.time()
print("DR Production amount (CVXPY + cvxRiskOpt + CVXPYgen): ", x.value)
print("Solve Time: %.3f ms" % (1000 * (t1 - t0)))

DR Production amount (CVXPY + cvxRiskOpt + CVXPYgen):  707.0193730343249
Solve Time: 0.199 ms

Notice that the result 707.019 matches that from earlier, but the solve time is significantly lower.