add coex model

README.md

@@ -45,6 +45,7 @@ List of models under `models` directory is listed below. Most models have a GAMS
 | [chenery](models/chenery)                                     | [chenery](https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_chenery.html)                                       | NLP           | Demo             |
 | [circuit](models/circuit)                                     | [circuit](https://www.gams.com/latest/noalib_ml/libhtml/noalib_circuit.html)                                         | NLP           | Demo             |
 | [clsp](models/clsp)                                           |                                                                                                                      | MIP           | Demo             |
+| [coex](models/coex)                                           | [coex](https://www.gams.com/latest/noalib_ml/libhtml/noalib_control3.html)                                           | NLP           | Demo             |
 | [control3](models/control3)                                   | [control3](https://www.gams.com/latest/noalib_ml/libhtml/noalib_control3.html)                                       | NLP           | Demo             |
 | [Corporate](models/Corporate)                                 | [Corporate](https://www.gams.com/latest/finlib_ml/libhtml/finlib_Corporate.html)                                     | LP            | Requires license |
 | [cpa](models/cpa)                                             | [cpa](https://www.gams.com/latest/noalib_ml/libhtml/noalib_cpa.html)                                                 | NLP           | Demo             |

models/coex/coex.py

@@ -0,0 +1,121 @@
+"""
+## GAMSSOURCE: https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_coex.html
+## LICENSETYPE: Demo
+## MODELTYPE: MIP
+## KEYWORDS: mixed integer linear programming, mathematical games, combinatorial optimization, peaceably coexisting armies of queens
+
+Peacefully Coexisting Armies of Queens (COEX)
+
+Two armies of queens (black and white) peacefully coexist on a
+chessboard when they are placed on the board in such a way that
+no two queens from opposing armies can attack each other. The
+problem is to find the maximum two equal-sized armies.
+
+Bosch, R, Mind Sharpener. OPTIMA MPS Newsletter (2000).
+"""
+
+from sys import stdout
+
+import gamspy.math as gpm
+from gamspy import (
+    Alias,
+    Container,
+    Equation,
+    Model,
+    Ord,
+    Problem,
+    Sense,
+    Set,
+    Sum,
+    Variable,
+    VariableType,
+)
+
+m = Container()
+
+i = Set(
+    m,
+    "i",
+    description="size of chess board",
+    records=[str(i + 1) for i in range(8)],
+)
+j, ii, jj = (Alias(m, ident, i) for ident in ["j", "ii", "jj"])
+
+M = Set(
+    m, "M", domain=[i, j, ii, jj], description="shared positions on the board"
+)
+M[i, j, ii, jj] = (
+    (Ord(i) == Ord(ii))
+    | (Ord(j) == Ord(jj))
+    | (gpm.abs(Ord(i) - Ord(ii)) == gpm.abs(Ord(j) - Ord(jj)))
+)
+
+b = Variable(
+    m,
+    "b",
+    domain=[i, j],
+    type=VariableType.BINARY,
+    description="square occupied by a black queen",
+)
+w = Variable(
+    m,
+    "w",
+    domain=[i, j],
+    type=VariableType.BINARY,
+    description="square occupied by a white queen",
+)
+
+tot = Variable(m, "tot", description="total queens in each army")
+
+eq1 = Equation(
+    m, "eq1", domain=[i, j, ii, jj], description="keep armies at peace"
+)
+eq1[M[i, j, ii, jj]] = b[i, j] + w[ii, jj] <= 1
+eq2 = Equation(
+    m,
+    "eq2",
+    description="add up all the black queens",
+    definition=tot == Sum((i, j), b[i, j]),
+)
+eq3 = Equation(
+    m,
+    "eq3",
+    description="add up all the white queens",
+    definition=tot == Sum((i, j), w[i, j]),
+)
+
+armies = Model(
+    m,
+    "armies",
+    problem=Problem.MIP,
+    equations=m.getEquations(),
+    sense=Sense.MAX,
+    objective=tot,
+)
+armies.solve(output=stdout)
+# Display solution
+print(f"Army size: {int(armies.objective_value)}")
+
+
+def queen_pos(var: Variable):
+    res_dict = var.toDict()
+    return [
+        (int(coord[0]), int(coord[1]))
+        for coord, v in res_dict.items()
+        if v > 0
+    ]
+
+
+bpos, wpos = queen_pos(b), queen_pos(w)
+print(f"Black queen positions: {bpos}\nWhite queen positions: {wpos}")
+
+
+# Validate solution
+def can_attack(ar, ac, br, bc):
+    return ar == br or ac == bc or abs(ar - br) == abs(ac - bc)
+
+
+assert int(armies.objective_value) == 9
+for ar, ac in bpos:
+    for br, bc in wpos:
+        assert not can_attack(ar, ac, br, bc)

models/refrigeration/refrigeration.py

@@ -18,142 +18,86 @@
 
 from __future__ import annotations
 
-from gamspy import Container, Equation, Model, Parameter, Variable
+from gamspy import Container, Equation, Model, Parameter, Sense, Variable
 
 
 def main():
     m = Container()
 
     # VARIABLES #
-    x1 = Variable(m, name="x1")
-    x2 = Variable(m, name="x2")
-    x3 = Variable(m, name="x3")
-    x4 = Variable(m, name="x4")
-    x5 = Variable(m, name="x5")
-    x6 = Variable(m, name="x6")
-    x7 = Variable(m, name="x7")
-    x8 = Variable(m, name="x8")
-    x9 = Variable(m, name="x9")
-    x10 = Variable(m, name="x10")
-    x11 = Variable(m, name="x11")
-    x12 = Variable(m, name="x12")
-    x13 = Variable(m, name="x13")
-    x14 = Variable(m, name="x14")
+    x = [Variable(m, name=f"x{i+1}") for i in range(14)]
 
     # EQUATIONS #
-    e1 = Equation(m, name="e1", type="regular")
-    e2 = Equation(m, name="e2", type="regular")
-    e3 = Equation(m, name="e3", type="regular")
-    e4 = Equation(m, name="e4", type="regular")
-    e5 = Equation(m, name="e5", type="regular")
-    e6 = Equation(m, name="e6", type="regular")
-    e7 = Equation(m, name="e7", type="regular")
-    e8 = Equation(m, name="e8", type="regular")
-    e9 = Equation(m, name="e9", type="regular")
-    e10 = Equation(m, name="e10", type="regular")
-    e11 = Equation(m, name="e11", type="regular")
-    e12 = Equation(m, name="e12", type="regular")
-    e13 = Equation(m, name="e13", type="regular")
-    e14 = Equation(m, name="e14", type="regular")
-    e15 = Equation(m, name="e15", type="regular")
+    e = [Equation(m, name=f"e{i+1}") for i in range(15)]
 
     # Objective function to be minimized:
     eobj = (
-        63098.88 * x2 * x4 * x12
-        + 5441.5 * x12 * x2**2
-        + 115055.5 * x6 * (x2**1.664)
-        + 6172.27 * x6 * x2**2
-        + 63098.88 * x1 * x3 * x11
-        + 5441.5 * x11 * x1**2
-        + 115055.5 * x5 * (x1**1.664)
-        + 6172.27 * x5 * x1**2
-        + 140.53 * x1 * x11
-        + 281.29 * x3 * x11
-        + 70.26 * x1**2
-        + 281.29 * x1 * x3
-        + 281.29 * x3**2
-        + 14437 * (x8**1.8812) * (x12**0.3424) * x7 * x10 * (x1**2) / (x9 * x14)
-        + 20470.2 * (x7**2.893) * (x11 * 0.316) * (x1 * 82)
+        63098.88 * x[1] * x[3] * x[11]
+        + 5441.5 * x[11] * x[1] ** 2
+        + 115055.5 * x[5] * (x[1] ** 1.664)
+        + 6172.27 * x[5] * x[1] ** 2
+        + 63098.88 * x[0] * x[2] * x[10]
+        + 5441.5 * x[10] * x[0] ** 2
+        + 115055.5 * x[4] * (x[0] ** 1.664)
+        + 6172.27 * x[4] * x[0] ** 2
+        + 140.53 * x[0] * x[10]
+        + 281.29 * x[2] * x[10]
+        + 70.26 * x[0] ** 2
+        + 281.29 * x[0] * x[2]
+        + 281.29 * x[2] ** 2
+        + 14437
+        * (x[7] ** 1.8812)
+        * (x[11] ** 0.3424)
+        * x[6]
+        * x[9]
+        * (x[0] ** 2)
+        / (x[8] * x[13])
+        + 20470.2 * (x[6] ** 2.893) * (x[10] * 0.316) * (x[0] * 82)
     )
 
     # Constaints:
-    e1[...] = 1.524 / x7 <= 1
-    e2[...] = 1.524 / x8 <= 1
-    e3[...] = 0.07789 * x1 - 2 * x9 / x7 <= 1
-    e4[...] = 7.05305 * (x1**2) * x10 / (x2 * x8 * x9 * x14) <= 1
-    e5[...] = 0.0833 * x14 / x13 <= 1
-    e6[...] = (
-        47.136 * x12 * (x2**0.333) / x10
-        - 1.333 * x8 * (x13**2.1195)
-        + 62.08 * (x13**2.1195) * (x8**0.2) / (x10 * x12)
+    e[0][...] = 1.524 / x[6] <= 1
+    e[1][...] = 1.524 / x[7] <= 1
+    e[2][...] = 0.07789 * x[0] - 2 * x[8] / x[6] <= 1
+    e[3][...] = (
+        7.05305 * (x[0] ** 2) * x[9] / (x[1] * x[7] * x[8] * x[13]) <= 1
+    )
+    e[4][...] = 0.0833 * x[13] / x[12] <= 1
+    e[5][...] = (
+        47.136 * x[11] * (x[1] ** 0.333) / x[9]
+        - 1.333 * x[7] * (x[12] ** 2.1195)
+        + 62.08 * (x[12] ** 2.1195) * (x[7] ** 0.2) / (x[9] * x[11])
         <= 1
     )
-    e7[...] = 0.04771 * x10 * (x8**1.8812) * (x12**0.3424) <= 1
-    e8[...] = 0.0488 * x9 * (x7**1.893) * (x11**0.316) <= 1
-    e9[...] = 0.0099 * x1 / x3 <= 1
-    e10[...] = 0.0193 * x2 / x4 <= 1
-    e11[...] = 0.0298 * x1 / x5 <= 1
-    e12[...] = 0.056 * x2 / x6 <= 1
-    e13[...] = 2 / x9 <= 1
-    e14[...] = 2 / x10 <= 1
-    e15[...] = x12 / x11 <= 1
+    e[6][...] = 0.04771 * x[9] * (x[7] ** 1.8812) * (x[11] ** 0.3424) <= 1
+    e[7][...] = 0.0488 * x[8] * (x[6] ** 1.893) * (x[10] ** 0.316) <= 1
+    e[8][...] = 0.0099 * x[0] / x[2] <= 1
+    e[9][...] = 0.0193 * x[1] / x[3] <= 1
+    e[10][...] = 0.0298 * x[0] / x[4] <= 1
+    e[11][...] = 0.056 * x[1] / x[5] <= 1
+    e[12][...] = 2 / x[8] <= 1
+    e[13][...] = 2 / x[9] <= 1
+    e[14][...] = x[11] / x[10] <= 1
 
     # Bounds on variables:
-    x1.lo[...] = 0.001
-    x1.up[...] = 5
-    x2.lo[...] = 0.001
-    x2.up[...] = 5
-    x3.lo[...] = 0.001
-    x3.up[...] = 5
-    x4.lo[...] = 0.001
-    x4.up[...] = 5
-    x5.lo[...] = 0.001
-    x5.up[...] = 5
-    x6.lo[...] = 0.001
-    x6.up[...] = 5
-    x7.lo[...] = 0.001
-    x7.up[...] = 5
-    x8.lo[...] = 0.001
-    x8.up[...] = 5
-    x9.lo[...] = 0.001
-    x9.up[...] = 5
-    x10.lo[...] = 0.001
-    x10.up[...] = 5
-    x11.lo[...] = 0.001
-    x11.up[...] = 5
-    x12.lo[...] = 0.001
-    x12.up[...] = 5
-    x13.lo[...] = 0.001
-    x13.up[...] = 5
-    x14.lo[...] = 0.001
-    x14.up[...] = 5
+    for v in m.getVariables():
+        v.lo = 0.001
+        v.up = 5
 
     refrigeration = Model(
         m,
         name="refrigeration",
         equations=m.getEquations(),
         problem="nlp",
-        sense="MIN",
+        sense=Sense.MIN,
         objective=eobj,
     )
     refrigeration.solve()
 
     # REPORTING PARAMETER
     rep = Parameter(m, name="rep", domain=["*", "*"])
-    rep["x1", "value"] = x1.l
-    rep["x2", "value"] = x2.l
-    rep["x3", "value"] = x3.l
-    rep["x4", "value"] = x4.l
-    rep["x5", "value"] = x5.l
-    rep["x6", "value"] = x6.l
-    rep["x7", "value"] = x7.l
-    rep["x8", "value"] = x8.l
-    rep["x9", "value"] = x9.l
-    rep["x10", "value"] = x10.l
-    rep["x11", "value"] = x11.l
-    rep["x12", "value"] = x12.l
-    rep["x13", "value"] = x13.l
-    rep["x14", "value"] = x14.l
+    for i in range(14):
+        rep[f"x{i+1}", "value"] = x[i].l
 
     print(
         "Objective Function Value: ",