Robust linear programming as an SOCP (ROBUSTLP)#

robustlp.py

"""
## 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
from gamspy import Container
from gamspy import Equation
from gamspy import Model
from gamspy import Parameter
from gamspy import Problem
from gamspy import Sense
from gamspy import Set
from gamspy import Sum
from gamspy import Variable
from gamspy.math import uniform


def main():
    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
        delayed_execution=int(os.getenv("DELAYED_EXECUTION", False)),
    )

    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()