8.1 Robust BIM microchip production problem#
Preamble: Install Pyomo and solvers#
HiGHS solves the linear and mixed-integer linear models. SCIP solves both versions of the final algebraic model with ball uncertainty, preserving integrality when requested. Both are open source. The Colab installation below supplies their Python interfaces; locally, install them with %pip install "pyomo>=6.10.1" highspy pyscipopt.
The optional comparison using Pyomo’s kernel conic components requires a separately installed and licensed MOSEK. Leave RUN_MOSEK_COMPARISON false to run the open-source route only. License files and credentials are private and are not part of this notebook.
import sys
if "google.colab" in sys.modules:
%pip install -q "pyomo>=6.10.1" highspy pyscipopt
import pyomo.environ as pyo
solver = "appsi_highs"
NLO_solver = "scip_direct"
MINLO_solver = "scip_direct"
SOLVER = pyo.SolverFactory(solver)
NLO_SOLVER = pyo.SolverFactory(NLO_solver)
MINLO_SOLVER = pyo.SolverFactory(MINLO_solver)
assert SOLVER.available(), f"Solver {solver} is not available."
assert NLO_SOLVER.available(), f"Solver {NLO_solver} is not available."
assert MINLO_SOLVER.available(), f"Solver {MINLO_solver} is not available."
RUN_MOSEK_COMPARISON = False
if RUN_MOSEK_COMPARISON:
CONIC_SOLVER = pyo.SolverFactory("mosek_direct")
assert CONIC_SOLVER.available(), "Install MOSEK and configure its license first."
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.patches as patches
import fractions
Original BIM production planning model#
The full description of the BIM production problem, can be found here. The resulting linear optimization problem was formulated as follows:
chips = ["logic", "memory"]
profits = {"logic": 12, "memory": 9}
copper = {"logic": 4, "memory": 2}
m = pyo.ConcreteModel("BIM basic problem")
m.chips = pyo.Set(initialize=chips)
m.x = pyo.Var(m.chips, within=pyo.NonNegativeReals)
m.profit = pyo.Objective(
expr=pyo.quicksum([profits[c] * m.x[c] for c in m.chips]),
sense=pyo.maximize,
)
m.silicon = pyo.Constraint(expr=m.x["logic"] <= 1000)
m.gemanium = pyo.Constraint(expr=m.x["memory"] <= 1500)
m.plastic = pyo.Constraint(expr=pyo.quicksum([m.x[c] for c in m.chips]) <= 1750)
m.copper = pyo.Constraint(
expr=pyo.quicksum(copper[c] * m.x[c] for c in m.chips) <= 4800
)
solve_results = SOLVER.solve(m)
pyo.assert_optimal_termination(solve_results)
print(
f"The optimal solution is x={[round(pyo.value(m.x[c]),3) for c in m.chips]} and yields a profit of {pyo.value(m.profit):.2f}"
)
The optimal solution is x=[650.0, 1100.0] and yields a profit of 17700.00
def ShowDuals(model):
print("The dual variable corresponding to:")
for c in model.component_objects(pyo.Constraint, active=True):
print(
f"- the constraint on {c} is equal to {str(fractions.Fraction(model.dual[c]))}"
)
m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT)
solve_results = SOLVER.solve(m)
pyo.assert_optimal_termination(solve_results)
ShowDuals(m)
The dual variable corresponding to:
- the constraint on silicon is equal to 0
- the constraint on gemanium is equal to 0
- the constraint on plastic is equal to 6
- the constraint on copper is equal to 3/2
Robust BIM production planning models#
Suppose now that there is uncertainty affecting the microchip production at BIM. Specifically, the company notices that not the amount of copper needed for the two types of microchips is not exactly 4 and 2 gr, but varies due to some external factors affecting the production process. How does this uncertainty affect the optimal production plan?
To get a feeling for what happens, let us first perform some simulations and data analysis on them. We start by simulating a sample of \(n=2000\) observed copper consumption pairs for the production of f logic chips and g memory chips. The amounts vary around the original values, 4 gr and 2 gr, respectively, according to two independent lognormal distributions.
plt.rcParams.update({"font.size": 12})
seed = 0
rng = np.random.default_rng(seed)
n = 2000
f = rng.lognormal(np.log(4.0), 0.005, n)
g = rng.lognormal(np.log(2.0), 0.005, n)
plt.figure()
plt.plot(f, g, ".")
plt.xlabel("Copper gr needed for logic chips")
plt.ylabel("Copper gr needed for memory chips")
plt.show()
Box uncertainty for copper consumption#
A very simple and somehow naive uncertainty set can be the minimal box that contains all the simulated data.
plt.figure()
plt.plot(f, g, ".")
currentAxis = plt.gca()
currentAxis.add_patch(
patches.Rectangle(
(min(f), min(g)),
max(f) - min(f),
max(g) - min(g),
fill=False,
color="r",
)
)
plt.xlabel("Grams of copper needed for logic chips")
plt.ylabel("Grams of copper needed for memory chips")
plt.tight_layout()
plt.show()
# calculate the upper and lower bounds for each uncertain parameter
lower = {"logic": min(f), "memory": min(g)}
upper = {"logic": max(f), "memory": max(g)}
print("Lower bounds", lower)
print("Upper bounds", upper)
Lower bounds {'logic': np.float64(3.922766922829344), 'memory': np.float64(1.9701110863753781)}
Upper bounds {'logic': np.float64(4.061793174956137), 'memory': np.float64(2.0328386701386703)}
Using this empirical box uncertainty set, we can consider the following robust variant of their optimization model:
The above model has an infinite number of constraints, one for every realization of the uncertain coefficients \(z\). However, using linear duality, we can deal with this and obtain a robustified linear optimization problem that we can solve.
Robust counterpart of box uncertainty#
The first thing to notice is that the copper consumption is modeled by constraints that are equivalent to bounding the following optimization problem:
or
Now we use linear duality to realize that the above is equivalent to:
and the constraint imposed by the last problem is equivalent to:
The only thing we need to do is add the new auxiliary variables and constraints to the original model and implement them in Pyomo.
def BIMWithBoxUncertainty(lower, upper, domain=pyo.NonNegativeReals):
m = pyo.ConcreteModel("BIM with Box Uncertainty")
m.chips = pyo.Set(initialize=chips)
m.x = pyo.Var(m.chips, within=domain)
m.profit = pyo.Objective(
expr=sum([profits[c] * m.x[c] for c in m.chips]), sense=pyo.maximize
)
m.silicon = pyo.Constraint(expr=m.x["logic"] <= 1000)
m.germanium = pyo.Constraint(expr=m.x["memory"] <= 1500)
m.plastic = pyo.Constraint(expr=sum([m.x[c] for c in m.chips]) <= 1750)
m.y = pyo.Var(m.chips, domain=pyo.NonNegativeReals)
m.w = pyo.Var(m.chips, domain=pyo.NonNegativeReals)
m.robustcopper = pyo.Constraint(
expr=sum([upper[c] * m.y[c] - lower[c] * m.w[c] for c in m.chips]) <= 4800
)
@m.Constraint(m.chips)
def PerVariable(m, c):
return m.x[c] == m.y[c] - m.w[c]
return m
m = BIMWithBoxUncertainty(lower, upper)
solve_results = SOLVER.solve(m)
pyo.assert_optimal_termination(solve_results)
print(
f"The optimal solution is x={[round(pyo.value(m.x[c]),3) for c in m.chips]} and yields a profit of {pyo.value(m.profit):.2f}"
)
The optimal solution is x=[612.4, 1137.6] and yields a profit of 17587.20
We may want to impose the box uncertainty set to be symmetric with respect to the nominal values and just choose its width \(\delta\). This leads to a different optimal robust solution.
# The parameter delta allows you to tune the amount of uncertainty.
# In particular, if you take delta=0, you obtain the same result as the nominal model.
delta = 0.05
def BIMWithSymmetricalBoxUncertainty(delta, domain=pyo.NonNegativeReals):
lower = {chip: copper[chip] - delta for chip in chips}
upper = {chip: copper[chip] + delta for chip in chips}
return BIMWithBoxUncertainty(lower, upper, domain=domain)
m = BIMWithSymmetricalBoxUncertainty(delta)
solve_results = SOLVER.solve(m)
pyo.assert_optimal_termination(solve_results)
print(
f"The optimal solution is x={[round(pyo.value(m.x[c]),3) for c in m.chips]} and yields a profit of {pyo.value(m.profit):.2f}"
)
The optimal solution is x=[606.25, 1143.75] and yields a profit of 17568.75
Integer solution variant#
The original BIM model gave integer solutions, but not the robust version. If we need integer solutions then we should impose that to the nature of the variables, which in this case of box uncertainty is easy to do since the model remains linear, although it will be mixed integer.
m = BIMWithBoxUncertainty(lower, upper, domain=pyo.NonNegativeIntegers)
solve_results = SOLVER.solve(m)
pyo.assert_optimal_termination(solve_results)
print(
f"The optimal solution is x={[round(pyo.value(m.x[c]),3) for c in m.chips]} and yields a profit of {pyo.value(m.profit):.2f}"
)
The optimal solution is x=[612.0, 1138.0] and yields a profit of 17586.00
Let us see how the optimal solution behave as we vary the width of the box uncertainty set \(\delta\) from 0 to 0.5.
df = pd.DataFrame()
for delta in np.linspace(0, 0.5, 21):
m = BIMWithSymmetricalBoxUncertainty(delta, domain=pyo.NonNegativeIntegers)
solve_results = SOLVER.solve(m)
pyo.assert_optimal_termination(solve_results)
results = [pyo.value(m.profit)] + [pyo.value(m.x[i]) for i in m.chips]
df.at[delta, "profit"] = results[0]
df.at[delta, chips[0]] = results[1]
df.at[delta, chips[1]] = results[2]
df
| profit | logic | memory | |
|---|---|---|---|
| 0.000 | 17700.0 | 650.0 | 1100.0 |
| 0.025 | 17634.0 | 628.0 | 1122.0 |
| 0.050 | 17568.0 | 606.0 | 1144.0 |
| 0.075 | 17502.0 | 584.0 | 1166.0 |
| 0.100 | 17436.0 | 562.0 | 1188.0 |
| 0.125 | 17370.0 | 540.0 | 1210.0 |
| 0.150 | 17304.0 | 518.0 | 1232.0 |
| 0.175 | 17238.0 | 496.0 | 1254.0 |
| 0.200 | 17175.0 | 475.0 | 1275.0 |
| 0.225 | 17109.0 | 453.0 | 1297.0 |
| 0.250 | 17043.0 | 431.0 | 1319.0 |
| 0.275 | 16977.0 | 409.0 | 1341.0 |
| 0.300 | 16911.0 | 387.0 | 1363.0 |
| 0.325 | 16845.0 | 365.0 | 1385.0 |
| 0.350 | 16779.0 | 343.0 | 1407.0 |
| 0.375 | 16713.0 | 321.0 | 1429.0 |
| 0.400 | 16650.0 | 300.0 | 1450.0 |
| 0.425 | 16584.0 | 278.0 | 1472.0 |
| 0.450 | 16518.0 | 256.0 | 1494.0 |
| 0.475 | 16416.0 | 243.0 | 1500.0 |
| 0.500 | 16296.0 | 233.0 | 1500.0 |
We can visualize how these quantities change as a function of \(\delta\):
df[["profit"]].plot()
plt.ylim([16001, 17999])
plt.xlabel(r"Margin $\delta$ of the uncertainty box")
plt.ylabel("Profit")
plt.show()
df[["logic", "memory"]].plot()
plt.xlabel(r"Margin $\delta$ of the uncertainty box")
plt.ylabel("Optimal number of produced chips")
plt.show()
Cardinality-constrained uncertainty set#
Let us now make different assumptions regarding the uncertainty related to the copper consumption. More specifically, we now assume that each uncertain coefficient \(z_j\) may deviate by at most \(\pm \delta\) from the nominal value \(\bar{z}_j\) with a total budget \(\sum_j |y_j|\leq\Gamma\) of normalized deviations.
Here \(\delta\) is an absolute deviation in grams per chip, not a percentage of the nominal coefficient. For fractional \(\Gamma\), the budget can be shared fractionally; for example, \(\Gamma=1.5\) allows one full deviation and half of another. We use this same uncertainty set in the robust counterpart and the adversarial algorithm below.
Robust counterpart of cardinality-constrained uncertainty#
Lagrange duality yields the following modification to the problem as equivalent to the robust model stated above:
def BIMWithBudgetUncertainty(delta, gamma, domain=pyo.NonNegativeReals):
m = pyo.ConcreteModel("BIM with Budget Uncertainty")
m.chips = pyo.Set(initialize=chips)
m.x = pyo.Var(m.chips, domain=domain)
m.profit = pyo.Objective(
expr=sum([profits[c] * m.x[c] for c in m.chips]), sense=pyo.maximize
)
m.silicon = pyo.Constraint(expr=m.x["logic"] <= 1000)
m.germanium = pyo.Constraint(expr=m.x["memory"] <= 1500)
m.plastic = pyo.Constraint(expr=sum([m.x[c] for c in m.chips]) <= 1750)
m.t = pyo.Var(m.chips, domain=pyo.NonNegativeReals)
m.lam = pyo.Var(domain=pyo.NonNegativeReals)
m.robustcopper = pyo.Constraint(
expr=sum([copper[c] * m.x[c] for c in m.chips])
+ gamma * m.lam
+ sum(m.t[c] for c in m.chips)
<= 4800
)
@m.Constraint(m.chips)
def up_rule(m, c):
return m.t[c] >= delta * m.x[c] - m.lam
@m.Constraint(m.chips)
def down_rule(m, c):
return m.t[c] >= -delta * m.x[c] - m.lam
return m
m = BIMWithBudgetUncertainty(0.01, 2, domain=pyo.NonNegativeIntegers)
solve_results = SOLVER.solve(m)
pyo.assert_optimal_termination(solve_results)
print(
f"The optimal solution is x={[round(pyo.value(m.x[c]),3) for c in m.chips]} and yields a profit of {pyo.value(m.profit):.2f}"
)
The optimal solution is x=[641.0, 1109.0] and yields a profit of 17673.00
Adversarial approach for the budgeted uncertainty set#
Instead of adopting the approach of robust counterparts, we could also use the adversarial approach where we initially solve the problem for the nominal value of the data. Then, we iteratively search for scenarios that make the current solution violate the copper constraint, and re-solve the problem to take this scenario into account. To do so, we need to slightly modify our problem formulation function to allow for many scenarios for the parameter \(z\).
def BIMWithSetOfScenarios(
delta, Z=None, domain=pyo.NonNegativeReals
):
if Z is None:
Z = [{"logic": 0, "memory": 0}]
chips = ["logic", "memory"]
profits = {"logic": 12, "memory": 9}
copper = {"logic": 4, "memory": 2}
m = pyo.ConcreteModel("BIM basic problem")
m.chips = pyo.Set(initialize=chips)
m.scenarios = pyo.Set(initialize=range(len(Z)))
m.x = pyo.Var(m.chips, within=domain)
m.profit = pyo.Objective(
expr=pyo.quicksum([profits[c] * m.x[c] for c in m.chips]),
sense=pyo.maximize,
)
m.silicon = pyo.Constraint(expr=m.x["logic"] <= 1000)
m.gemanium = pyo.Constraint(expr=m.x["memory"] <= 1500)
m.plastic = pyo.Constraint(expr=pyo.quicksum([m.x[c] for c in m.chips]) <= 1750)
@m.Constraint(m.scenarios)
def balance(m, i):
z = Z[i]
return pyo.quicksum((copper[c] + delta * z[c]) * m.x[c] for c in m.chips) <= 4800
return m
We also need a function that for a given solution finds the worst-possible realization of the uncertainty restricted by the parameters \(\Gamma\) and \(\delta\). In other words, its role is to solve the following maximization problem for a given solution \((\bar{x}_1, \bar{x}_2)\):
Such a function is implemented below and takes as argument also the maximum magnitude of the individual deviations \(\delta\) and the total budget \(\Gamma\).
def BIMPessimization(x, delta, gamma):
chips = ["logic", "memory"]
copper = {"logic": 4, "memory": 2}
m = pyo.ConcreteModel("BIM pessimization problem")
m.chips = pyo.Set(initialize=chips)
m.z = pyo.Var(m.chips, within=pyo.Reals)
m.u = pyo.Var(m.chips, within=pyo.NonNegativeReals)
@m.Constraint(m.chips)
def absolute_value_1(m, i):
return m.z[i] <= m.u[i]
@m.Constraint(m.chips)
def absolute_value_2(m, i):
return -m.z[i] <= m.u[i]
@m.Constraint(m.chips)
def absolute_value_less_than_one(m, i):
return m.u[i] <= 1.0
m.budget_constraint = pyo.Constraint(
expr=pyo.quicksum([m.u[i] for i in m.chips]) <= gamma
)
m.violation = pyo.Objective(
expr=-4800
+ pyo.quicksum([(copper[c] + delta * m.z[c]) * x[c] for c in m.chips]),
sense=pyo.maximize,
)
return m
We wrap the two functions above into a loop of the adversarial approach, which begins with a non-perturbation assumption and gradually generates violating scenarios, reoptimizing until the maximum constraint violation is below a tolerable threshold.
# Parameters
adversarial_converged = False
stopping_precision = 1e-7
max_iterations = 20
adversarial_iterations = 0
delta = 0.2
gamma = 1.5
chips = ["logic", "memory"]
# Initialize the null scenario - no perturbation
Z = [{"logic": 0, "memory": 0}]
while (not adversarial_converged) and (adversarial_iterations < max_iterations):
# Building and solving the master problem
model = BIMWithSetOfScenarios(delta, Z, domain=pyo.NonNegativeIntegers)
solve_results = SOLVER.solve(model)
pyo.assert_optimal_termination(solve_results)
# Saving the current solution
x = {i: model.x[i]() for i in model.chips}
print(f"\nIteration #{adversarial_iterations}")
print(f"Current solution: ")
for c in chips:
print(f"x['{c}']= {x[c]:.2f}")
# Pessimization
m = BIMPessimization(x, delta, gamma)
solve_results = SOLVER.solve(m)
pyo.assert_optimal_termination(solve_results)
worst_z = {i: m.z[i]() for i in chips}
constraint_violation = m.violation()
# If pessimization yields no violation, stop the procedure, otherwise add a scenario and repeat
if constraint_violation < stopping_precision:
print("No violation found. Stopping the procedure.")
adversarial_converged = True
else:
print(
f"Violation found: z['logic'] = {worst_z['logic']}, z['memory'] = {worst_z['memory']}, "
f"constraint violation: {constraint_violation:6.2f}"
)
Z.append(worst_z)
adversarial_iterations += 1
if not adversarial_converged:
raise RuntimeError("The iteration limit was reached before robust feasibility was verified.")
counterpart = BIMWithBudgetUncertainty(delta, gamma, domain=pyo.NonNegativeIntegers)
counterpart_results = SOLVER.solve(counterpart)
pyo.assert_optimal_termination(counterpart_results)
print(f"Adversarial profit: {model.profit():.2f}")
print(f"Robust-counterpart profit: {counterpart.profit():.2f}")
assert abs(model.profit() - counterpart.profit()) < 1e-6
Iteration #0
Current solution:
x['logic']= 650.00
x['memory']= 1100.00
Violation found: z['logic'] = 0.5, z['memory'] = 1.0, constraint violation: 285.00
Iteration #1
Current solution:
x['logic']= 500.00
x['memory']= 1250.00
No violation found. Stopping the procedure.
Adversarial profit: 17250.00
Robust-counterpart profit: 17250.00
The last master solution is accepted only after pessimization confirms that its maximum copper-constraint violation is below the tolerance. Each master problem relaxes the full robust problem, so its optimal profit is an upper bound. Once its solution is also robustly feasible, those two facts certify optimality, up to the stated tolerances. The direct robust counterpart provides an additional check using exactly the same \(\delta\), \(\Gamma\) and integer domain.
For this two-dimensional budget set, an exact linear pessimization can return a vertex. There are finitely many vertices, so adding a new violated vertex each time gives finite convergence. A rounded output or an iteration limit alone is not a convergence certificate. For a ball uncertainty set there are infinitely many extreme points, which is one reason to prefer a conic reformulation when available.
We will now illustrate how to use conic optimization to solve the problem using robust counterparts for an ellipsoidal uncertainty set.
Ball uncertainty set#
Let us now make yet another different assumption regarding the uncertainty related to copper consumption. More specifically, we assume that the two uncertain coefficients \(z_1\) and \(z_2\) can vary in a 2-dimensional ball centered around the point \((\bar{z}_1,\bar{z}_2) = (4,2)\) and with radius \(r\).
Robust counterpart of ball uncertainty#
A straightforward reformulation leads to the equivalent constraint:
By defining \(y = 4800 - \bar{z}_1x_1 - \bar{z}_2x_2\) and \(w = r x\), we may write:
We can express the ball constraint using Pyomo’s kernel components, an alternative modeling interface with explicit conic constraints. The next function illustrates this representation. Its two solve cells are optional MOSEK comparisons; the subsequent pyomo.environ formulation supplies the open-source SCIP route for both continuous and integer production decisions.
import pyomo.kernel as pyk
def BIMWithBallUncertainty(radius, domain_type=pyk.RealSet):
idxChips = range(len(chips))
m = pyk.block()
m.x = pyk.variable_list()
for i in idxChips:
m.x.append(pyk.variable(lb=0, domain=domain_type))
m.profit = pyk.objective(
expr=sum(profits[chips[i]] * m.x[i] for i in idxChips),
sense=pyk.maximize,
)
m.silicon = pyk.constraint(expr=m.x[0] <= 1000)
m.germanium = pyk.constraint(expr=m.x[1] <= 1500)
m.plastic = pyk.constraint(expr=sum([m.x[i] for i in idxChips]) <= 1750)
m.y = pyk.variable(lb=0)
m.w = pyk.variable_list()
for i in idxChips:
m.w.append(pyk.variable(lb=0))
m.copper = pyk.constraint(
expr=m.y == 4800 - sum(copper[chips[i]] * m.x[i] for i in idxChips)
)
m.xtow = pyk.constraint_list()
for i in idxChips:
m.xtow.append(pyk.constraint(expr=m.w[i] == radius * m.x[i]))
from pyomo.core.kernel.conic import quadratic
m.robust = quadratic(m.y, m.w)
return m
The following optional comparison uses MOSEK to solve the continuous conic formulation. Enable RUN_MOSEK_COMPARISON after installing MOSEK and configuring its license.
radius = 0.05
if RUN_MOSEK_COMPARISON:
m = BIMWithBallUncertainty(radius)
results = CONIC_SOLVER.solve(m)
pyo.assert_optimal_termination(results)
print(
f"Solver: MOSEK, solver status:",
results.solver.status,
"and solver terminal condition:",
results.solver.termination_condition,
)
print(
f"The optimal solution is x={[round(pyk.value(m.x[i]), 3) for i in range(len(m.x))]} and yields a profit of {pyo.value(m.profit):.2f}"
)
else:
print("Optional continuous MOSEK comparison is disabled; SCIP results follow below.")
Optional continuous MOSEK comparison is disabled; SCIP results follow below.
With integer production variables, the conic model requires a solver that enforces integrality. This optional cell uses MOSEK; the algebraic model below uses SCIP.
if RUN_MOSEK_COMPARISON:
m = BIMWithBallUncertainty(radius, domain_type=pyk.IntegerSet)
results = CONIC_SOLVER.solve(m)
pyo.assert_optimal_termination(results)
print(
f"Solver: MOSEK, solver status:",
results.solver.status,
"and solver terminal condition:",
results.solver.termination_condition,
)
print(
f"The optimal solution is x={[round(pyk.value(m.x[i]), 3) for i in range(len(m.x))]} and yields a profit of {pyo.value(m.profit):.2f}"
)
else:
print("Optional integer MOSEK comparison is disabled; SCIP results follow below.")
Optional integer MOSEK comparison is disabled; SCIP results follow below.
Implementing second-order cones using pyomo.environ#
For \(t\geq0\), the inequality \(\|x\|_2\leq t\) is equivalent to \(\|x\|_2^2\leq t^2\). Although the squared expression is not a convex function of all its arguments, the feasible set with non-negative \(t\) is the second-order cone. This gives an algebraic model that SCIP can solve, with either continuous or integer production variables. Explicit variable bounds also help its global optimization algorithm.
def BIMWithBallUncertaintyAsSquaredSecondOrderCone(r, domain=pyo.NonNegativeReals):
m = pyo.ConcreteModel("BIM with Ball Uncertainty as SOC")
m.chips = pyo.Set(initialize=chips)
m.x = pyo.Var(m.chips, within=domain, bounds=lambda m, c: (0, 1000 if c == "logic" else 1500))
# the nonnegativity of this variable is essential!
m.y = pyo.Var(within=pyo.NonNegativeReals, bounds=(0, 4800))
m.profit = pyo.Objective(
expr=sum([profits[c] * m.x[c] for c in m.chips]), sense=pyo.maximize
)
m.silicon = pyo.Constraint(expr=m.x["logic"] <= 1000)
m.germanium = pyo.Constraint(expr=m.x["memory"] <= 1500)
m.plastic = pyo.Constraint(expr=sum([m.x[c] for c in m.chips]) <= 1750)
m.copper = pyo.Constraint(
expr=m.y == 4800 - sum(copper[c] * m.x[c] for c in m.chips)
)
m.robust = pyo.Constraint(expr=sum((r * m.x[c]) ** 2 for c in m.chips) <= m.y**2)
return m
for label, domain, chosen_solver in [
("Continuous", pyo.NonNegativeReals, NLO_SOLVER),
("Integer", pyo.NonNegativeIntegers, MINLO_SOLVER),
]:
m = BIMWithBallUncertaintyAsSquaredSecondOrderCone(radius, domain=domain)
results = chosen_solver.solve(m)
pyo.assert_optimal_termination(results)
values = [pyo.value(m.x[c]) for c in m.chips]
print(f"{label} SCIP solution: x={[round(v, 3) for v in values]}, profit = {m.profit():.2f}")
if domain is pyo.NonNegativeIntegers:
assert all(abs(v - round(v)) < 1e-6 for v in values)
Continuous SCIP solution: x=[617.755, 1132.245], profit = 17603.26
Integer SCIP solution: x=[617.0, 1133.0], profit = 17601.00
m = BIMWithBallUncertaintyAsSquaredSecondOrderCone(
radius, domain=pyo.NonNegativeIntegers
)
results = MINLO_SOLVER.solve(m)
print(
f"Solver: {MINLO_solver}, solver status:",
results.solver.status,
"and solver terminal condition:",
results.solver.termination_condition,
)
print(
f"The optimal solution is x={[round(pyo.value(m.x[c]),3) for c in m.chips]} and yields a profit of {pyo.value(m.profit):.2f}"
)
Solver: scip_direct, solver status: ok and solver terminal condition: optimal
The optimal solution is x=[617.0, 1133.0] and yields a profit of 17601.00