"""
## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_robustlp.html
## LICENSETYPE: Demo
## MODELTYPE: LP, QCP
## KEYWORDS: linear programming, quadratic constraint programming, robust optimization, second order cone programming
Robust linear programming as an SOCP (ROBUSTLP)
Consider a linear optimization problem of the form
min_x c^Tx s.t. a_i^Tx <= b_i, i=1,..,m.
In practice, the coefficient vectors a_i may not be known perfectly,
as they are subject to noise. Assume that we only know that a_i in E_i,
where E_i are given ellipsoids. In robust optimization, we seek to minimize
the original objective, but we insist that each constraint be satisfied,
irrespective of the choice of the corresponding vector a_i in E_i.
We obtain the second-order cone optimization problem
min_x c^Tx s.t. a'_i^Tx + ||R_i^Tx|| <= b_i, i=1,..,m,
where E_i = { a'_i + R_iu | ||u|| <= 1}. In the above, we observe that
the feasible set is smaller than the original one, due to the terms involving
the l_2-norms.
The figure above illustrates the kind of feasible set one obtains in a
particular instance of the above problem, with spherical uncertainties
(that is, all the ellipsoids are spheres, R_i = rho I for some rho >0).
We observe that the robust feasible set is indeed contained in the original
polyhedron.
In this particular example we allow coefficients A(i,*) to vary in an
ellipsoid. The robust LP is reformulated as a SOCP.
Contributed by Michael Ferris, University of Wisconsin, Madison
Lobo, M S, Vandenberghe, L, Boyd, S, and Lebret, H, Applications of
Second Order Cone Programming. Linear Algebra and its Applications,
Special Issue on Linear Algebra in Control, Signals and Image
Processing. 284 (November, 1998).
"""
from __future__ import annotations
import os
from gamspy import (
Alias,
Container,
Equation,
Model,
Parameter,
Problem,
Sense,
Set,
Sum,
Variable,
)
from gamspy.math import uniform
def main():
m = Container(
system_directory=os.getenv("SYSTEM_DIRECTORY", None),
)
mu = 1e-2
# Set
i = Set(m, name="i", records=[f"{idx}" for idx in range(1, 8)])
j = Set(m, name="j", records=[f"{idx}" for idx in range(1, 5)])
# Data
a = Parameter(m, name="a", domain=[i, j])
b = Parameter(m, name="b", domain=i)
c = Parameter(m, name="c", domain=j)
b[i] = 1
c[j] = -1
a[i, j] = uniform(0, 1)
# Variable
x = Variable(m, name="x", domain=j)
# Equation
cons = Equation(m, name="cons", domain=i)
defobj = Sum(j, c[j] * x[j])
cons[i] = Sum(j, a[i, j] * x[j]) <= b[i]
lpmod = Model(
m,
name="lpmod",
equations=[cons],
problem="LP",
sense=Sense.MIN,
objective=defobj,
)
lpmod.solve()
results = Parameter(m, name="results", domain=["*", "*"])
results["lp", j] = x.l[j]
results["lp", "obj"] = lpmod.objective_value
lmbda = Variable(m, name="lambda", type="positive", domain=j)
gamma = Variable(m, name="gamma", type="positive", domain=j)
lpcons = Equation(m, name="lpcons", domain=i)
defdual = Equation(m, name="defdual", domain=j)
lpcons[i] = (
mu * Sum(j, lmbda[j] + gamma[j]) + Sum(j, a[i, j] * x[j]) <= b[i]
)
defdual[j] = lmbda[j] - gamma[j] == x[j]
lproblp = Model(
m,
name="lproblp",
equations=[lpcons, defdual],
problem="LP",
sense=Sense.MIN,
objective=defobj,
)
lproblp.solve()
results["roblp", j] = x.l[j]
results["roblp", "obj"] = lproblp.objective_value
k = Alias(m, name="k", alias_with=j)
p = Parameter(m, name="p", domain=[i, j, k])
p[i, j, j] = mu
y = Variable(m, name="y", domain=i)
v = Variable(m, name="v", domain=[i, k])
defrhs = Equation(m, name="defrhs", domain=i)
defv = Equation(m, name="defv", domain=[i, k])
socpqcpcons = Equation(m, name="socpqcpcons", domain=i)
defrhs[i] = y[i] == b[i] - Sum(j, a[i, j] * x[j])
defv[i, k] = v[i, k] == Sum(j, p[i, j, k] * x[j])
socpqcpcons[i] = y[i] ** 2 >= Sum(k, v[i, k] ** 2)
roblpqcp = Model(
m,
name="roblpqcp",
equations=[socpqcpcons, defrhs, defv],
problem=Problem.QCP,
sense=Sense.MIN,
objective=defobj,
)
y.lo[i] = 0
roblpqcp.solve()
results["qcp", j] = x.l[j]
results["qcp", "obj"] = roblpqcp.objective_value
print(results.records)
if __name__ == "__main__":
main()
