Skip to content

Commit 3886f4a

Browse files
Joao-DionisioJoao-Dionisio
committed
Clean up code and example
1 parent 06b4e72 commit 3886f4a

File tree

2 files changed

+42
-142
lines changed

2 files changed

+42
-142
lines changed

‎examples/finished/plot_primal_dual_evolution.py‎

Lines changed: 36 additions & 122 deletions
Original file line numberDiff line numberDiff line change
@@ -1,140 +1,54 @@
1+
"""
2+
This example show how to retrieve the primal and dual solutions during the optimization process
3+
and plot them as a function of time. The model is about gas transportation and can be found in
4+
PySCIPOpt/tests/helpers/utils.py
5+
6+
It makes use of the attach_primal_dual_evolution_eventhdlr recipe.
7+
8+
Requires matplotlib, and may require PyQt6 to show the plot.
9+
"""
10+
111
from pyscipopt import Model
212

313
def plot_primal_dual_evolution(model: Model):
414
try:
515
from matplotlib import pyplot as plt
616
except ImportError:
7-
raise("matplotlib is required to plot the solution. Try running `pip install matplotlib` in the command line.")
8-
9-
time_primal, val_primal = zip(*model.data["primal_log"])
17+
raise ImportError("matplotlib is required to plot the solution. Try running `pip install matplotlib` in the command line.\
18+
You may also need to install PyQt6 to show the plot.")
19+
20+
assert model.data["primal_log"], "Could not find any feasible solutions"
21+
time_primal, val_primal = map(list,zip(*model.data["primal_log"]))
22+
time_dual, val_dual = map(list,zip(*model.data["dual_log"]))
23+
24+
25+
if time_primal[-1] < time_dual[-1]:
26+
time_primal.append(time_dual[-1])
27+
val_primal.append(val_primal[-1])
28+
29+
if time_primal[-1] > time_dual[-1]:
30+
time_dual.append(time_primal[-1])
31+
val_dual.append(val_dual[-1])
32+
1033
plt.plot(time_primal, val_primal, label="Primal bound")
11-
time_dual, val_dual = zip(*model.data["dual_log"])
1234
plt.plot(time_dual, val_dual, label="Dual bound")
1335

1436
plt.legend(loc="best")
1537
plt.show()
1638

17-
18-
def gastrans_model():
19-
GASTEMP = 281.15
20-
RUGOSITY = 0.05
21-
DENSITY = 0.616
22-
COMPRESSIBILITY = 0.8
23-
nodes = [
24-
# name supplylo supplyup pressurelo pressureup cost
25-
("Anderlues", 0.0, 1.2, 0.0, 66.2, 0.0), # 0
26-
("Antwerpen", None, -4.034, 30.0, 80.0, 0.0), # 1
27-
("Arlon", None, -0.222, 0.0, 66.2, 0.0), # 2
28-
("Berneau", 0.0, 0.0, 0.0, 66.2, 0.0), # 3
29-
("Blaregnies", None, -15.616, 50.0, 66.2, 0.0), # 4
30-
("Brugge", None, -3.918, 30.0, 80.0, 0.0), # 5
31-
("Dudzele", 0.0, 8.4, 0.0, 77.0, 2.28), # 6
32-
("Gent", None, -5.256, 30.0, 80.0, 0.0), # 7
33-
("Liege", None, -6.385, 30.0, 66.2, 0.0), # 8
34-
("Loenhout", 0.0, 4.8, 0.0, 77.0, 2.28), # 9
35-
("Mons", None, -6.848, 0.0, 66.2, 0.0), # 10
36-
("Namur", None, -2.120, 0.0, 66.2, 0.0), # 11
37-
("Petange", None, -1.919, 25.0, 66.2, 0.0), # 12
38-
("Peronnes", 0.0, 0.96, 0.0, 66.2, 1.68), # 13
39-
("Sinsin", 0.0, 0.0, 0.0, 63.0, 0.0), # 14
40-
("Voeren", 20.344, 22.012, 50.0, 66.2, 1.68), # 15
41-
("Wanze", 0.0, 0.0, 0.0, 66.2, 0.0), # 16
42-
("Warnand", 0.0, 0.0, 0.0, 66.2, 0.0), # 17
43-
("Zeebrugge", 8.87, 11.594, 0.0, 77.0, 2.28), # 18
44-
("Zomergem", 0.0, 0.0, 0.0, 80.0, 0.0) # 19
45-
]
46-
arcs = [
47-
# node1 node2 diameter length active */
48-
(18, 6, 890.0, 4.0, False),
49-
(18, 6, 890.0, 4.0, False),
50-
(6, 5, 890.0, 6.0, False),
51-
(6, 5, 890.0, 6.0, False),
52-
(5, 19, 890.0, 26.0, False),
53-
(9, 1, 590.1, 43.0, False),
54-
(1, 7, 590.1, 29.0, False),
55-
(7, 19, 590.1, 19.0, False),
56-
(19, 13, 890.0, 55.0, False),
57-
(15, 3, 890.0, 5.0, True),
58-
(15, 3, 395.0, 5.0, True),
59-
(3, 8, 890.0, 20.0, False),
60-
(3, 8, 395.0, 20.0, False),
61-
(8, 17, 890.0, 25.0, False),
62-
(8, 17, 395.0, 25.0, False),
63-
(17, 11, 890.0, 42.0, False),
64-
(11, 0, 890.0, 40.0, False),
65-
(0, 13, 890.0, 5.0, False),
66-
(13, 10, 890.0, 10.0, False),
67-
(10, 4, 890.0, 25.0, False),
68-
(17, 16, 395.5, 10.5, False),
69-
(16, 14, 315.5, 26.0, True),
70-
(14, 2, 315.5, 98.0, False),
71-
(2, 12, 315.5, 6.0, False)
72-
]
73-
74-
model = Model()
75-
76-
# create flow variables
77-
flow = {}
78-
for arc in arcs:
79-
flow[arc] = model.addVar("flow_%s_%s" % (nodes[arc[0]][0], nodes[arc[1]][0]), # names of nodes in arc
80-
lb=0.0 if arc[4] else None) # no lower bound if not active
81-
82-
# pressure difference variables
83-
pressurediff = {}
84-
for arc in arcs:
85-
pressurediff[arc] = model.addVar("pressurediff_%s_%s" % (nodes[arc[0]][0], nodes[arc[1]][0]),
86-
# names of nodes in arc
87-
lb=None)
88-
89-
# supply variables
90-
supply = {}
91-
for node in nodes:
92-
supply[node] = model.addVar("supply_%s" % (node[0]), lb=node[1], ub=node[2], obj=node[5])
93-
94-
# square pressure variables
95-
pressure = {}
96-
for node in nodes:
97-
pressure[node] = model.addVar("pressure_%s" % (node[0]), lb=node[3] ** 2, ub=node[4] ** 2)
98-
99-
# node balance constrains, for each node i: outflows - inflows = supply
100-
for nid, node in enumerate(nodes):
101-
# find arcs that go or end at this node
102-
flowbalance = 0
103-
for arc in arcs:
104-
if arc[0] == nid: # arc is outgoing
105-
flowbalance += flow[arc]
106-
elif arc[1] == nid: # arc is incoming
107-
flowbalance -= flow[arc]
108-
else:
109-
continue
110-
111-
model.addCons(flowbalance == supply[node], name="flowbalance%s" % node[0])
112-
113-
# pressure difference constraints: pressurediff[node1 to node2] = pressure[node1] - pressure[node2]
114-
for arc in arcs:
115-
model.addCons(pressurediff[arc] == pressure[nodes[arc[0]]] - pressure[nodes[arc[1]]],
116-
"pressurediffcons_%s_%s" % (nodes[arc[0]][0], nodes[arc[1]][0]))
117-
118-
# pressure loss constraints:
119-
from math import log10
120-
for arc in arcs:
121-
coef = 96.074830e-15 * arc[2] ** 5 * (2.0 * log10(3.7 * arc[2] / RUGOSITY)) ** 2 / COMPRESSIBILITY / GASTEMP / \
122-
arc[3] / DENSITY
123-
if arc[4]: # active
124-
model.addCons(flow[arc] ** 2 + coef * pressurediff[arc] <= 0.0,
125-
"pressureloss_%s_%s" % (nodes[arc[0]][0], nodes[arc[1]][0]))
126-
else:
127-
model.addCons(flow[arc] * abs(flow[arc]) - coef * pressurediff[arc] == 0.0,
128-
"pressureloss_%s_%s" % (nodes[arc[0]][0], nodes[arc[1]][0]))
129-
130-
return model
131-
132-
13339
if __name__=="__main__":
13440
from pyscipopt.recipes.primal_dual_evolution import attach_primal_dual_evolution_eventhdlr
135-
41+
import os
42+
import sys
43+
44+
# just a way to import files from different folders, not important
45+
sys.path.append(os.path.abspath(os.path.join(os.path.dirname(__file__), '../../tests/helpers')))
46+
47+
from utils import gastrans_model
48+
13649
model = gastrans_model()
137-
model.attach_primal_dual_evolution_eventhdlr()
50+
model.data = {}
51+
attach_primal_dual_evolution_eventhdlr(model)
13852

13953
model.optimize()
14054
plot_primal_dual_evolution(model)

‎src/pyscipopt/recipes/primal_dual_evolution.py‎

Lines changed: 6 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -22,34 +22,20 @@ def eventinit(self): # we want to collect best primal solutions and best dual so
2222
def eventexec(self, event):
2323
# if a new best primal solution was found, we save when it was found and also its objective
2424
if event.getType() == SCIP_EVENTTYPE.BESTSOLFOUND:
25-
self.model.data["primal_log"].append((self.model.getSolvingTime(), self.model.getPrimalbound()))
25+
self.model.data["primal_log"].append([self.model.getSolvingTime(), self.model.getPrimalbound()])
2626

2727
if not self.model.data["dual_log"]:
28-
self.model.data["dual_log"].append((self.model.getSolvingTime(), self.model.getDualbound()))
28+
self.model.data["dual_log"].append([self.model.getSolvingTime(), self.model.getDualbound()])
2929

3030
if self.model.getObjectiveSense() == "minimize":
3131
if self.model.isGT(self.model.getDualbound(), self.model.data["dual_log"][-1][1]):
32-
self.model.data["dual_log"].append((self.model.getSolvingTime(), self.model.getDualbound()))
32+
self.model.data["dual_log"].append([self.model.getSolvingTime(), self.model.getDualbound()])
3333
else:
3434
if self.model.isLT(self.model.getDualbound(), self.model.data["dual_log"][-1][1]):
35-
self.model.data["dual_log"].append((self.model.getSolvingTime(), self.model.getDualbound()))
35+
self.model.data["dual_log"].append([self.model.getSolvingTime(), self.model.getDualbound()])
36+
3637

37-
def eventexitsol(self):
38-
if self.model.data["primal_log"][-1] and self.model.getPrimalbound() != self.model.data["primal_log"][-1][1]:
39-
self.model.data["primal_log"].append((self.model.getSolvingTime(), self.model.getPrimalbound()))
40-
41-
if not self.model.data["dual_log"]:
42-
self.model.data["dual_log"].append((self.model.getSolvingTime(), self.model.getDualbound()))
43-
44-
if self.model.getObjectiveSense() == "minimize":
45-
if self.model.isGT(self.model.getDualbound(), self.model.data["dual_log"][-1][1]):
46-
self.model.data["dual_log"].append((self.model.getSolvingTime(), self.model.getDualbound()))
47-
else:
48-
if self.model.isLT(self.model.getDualbound(), self.model.data["dual_log"][-1][1]):
49-
self.model.data["dual_log"].append((self.model.getSolvingTime(), self.model.getDualbound()))
50-
51-
52-
if not hasattr(model, "data"):
38+
if not hasattr(model, "data") or model.data==None:
5339
model.data = {}
5440

5541
model.data["primal_log"] = []

0 commit comments

Comments
 (0)