From 32e581ffefe1a89c2feb415dd104e97cc330950e Mon Sep 17 00:00:00 2001 From: mabdullahsoyturk Date: Fri, 16 Aug 2024 14:08:14 +0200 Subject: [PATCH] add coex model --- README.md | 1 + models/coex/coex.py | 121 ++++++++++++++++++++ models/refrigeration/refrigeration.py | 156 +++++++++----------------- 3 files changed, 172 insertions(+), 106 deletions(-) create mode 100644 models/coex/coex.py diff --git a/README.md b/README.md index 9ba1429..ff5b726 100644 --- a/README.md +++ b/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 | diff --git a/models/coex/coex.py b/models/coex/coex.py new file mode 100644 index 0000000..937310e --- /dev/null +++ b/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) diff --git a/models/refrigeration/refrigeration.py b/models/refrigeration/refrigeration.py index 23e93f8..3ac28fc 100644 --- a/models/refrigeration/refrigeration.py +++ b/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: ",