Scippy

SCIP

Solving Constraint Integer Programs

gastrans.c
Go to the documentation of this file.
1 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
2 /* */
3 /* This file is part of the program and library */
4 /* SCIP --- Solving Constraint Integer Programs */
5 /* */
6 /* Copyright (C) 2002-2022 Konrad-Zuse-Zentrum */
7 /* fuer Informationstechnik Berlin */
8 /* */
9 /* SCIP is distributed under the terms of the ZIB Academic License. */
10 /* */
11 /* You should have received a copy of the ZIB Academic License */
12 /* along with SCIP; see the file COPYING. If not visit scipopt.org. */
13 /* */
14 /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
15 
16 /**@file gastrans.c
17  * @brief Simple Gas Transportation Model
18  * @author Stefan Vigerske
19  *
20  * This example shows how to setup nonlinear constraints with signpower-expression when using SCIP as callable library.
21  * The example implements a model for the distribution of gas through a network of pipelines, which
22  * is formulated as a cost minimization subject to nonlinear flow-pressure relations, material balances,
23  * and pressure bounds. The Belgian gas network is used as an example.
24  *
25  * The model is taken from the GAMS model library:
26  * https://www.gams.com/latest/gamslib_ml/libhtml/gamslib_gastrans.html
27  *
28  * Original model source:
29  * @par
30  * D. de Wolf and Y. Smeers@n
31  * The Gas Transmission Problem Solved by and Extension of the Simplex Algorithm@n
32  * Management Science 46, 11 (2000), 1454-1465
33  */
34 
35 /*--+----1----+----2----+----3----+----4----+----5----+----6----+----7----+----8----+----9----+----0----+----1----+----2*/
36 
37 #include <stdio.h>
38 #include <math.h>
39 
40 #include "scip/scip.h"
41 #include "scip/scipdefplugins.h"
42 
43 /* Model parameters */
44 
45 /** node data structure */
46 typedef struct NodeData {
47  const char* name;
53 } NodeData;
54 
55 /** arc data structure */
56 typedef struct ArcData {
57  int node1;
58  int node2;
62 } ArcData;
63 
64 /** number of nodes (towns) */
65 #define nnodes 20
66 
67 /** number of arcs */
68 #define narcs 24
69 
70 /** value we use to represent infinity */
71 #define infinity 1e+20
72 
73 /** data of nodes */
74 static const NodeData nodedata[] =
75 {
76  /* name supplylo supplyup pressurelo pressureup cost */
77  {"Anderlues", 0.0, 1.2, 0.0, 66.2, 0.0 }, /* 0 */
78  {"Antwerpen", -infinity, -4.034, 30.0, 80.0, 0.0 }, /* 1 */
79  {"Arlon", -infinity, -0.222, 0.0, 66.2, 0.0 }, /* 2 */
80  {"Berneau", 0.0, 0.0, 0.0, 66.2, 0.0 }, /* 3 */
81  {"Blaregnies", -infinity, -15.616, 50.0, 66.2, 0.0 }, /* 4 */
82  {"Brugge", -infinity, -3.918, 30.0, 80.0, 0.0 }, /* 5 */
83  {"Dudzele", 0.0, 8.4, 0.0, 77.0, 2.28 }, /* 6 */
84  {"Gent", -infinity, -5.256, 30.0, 80.0, 0.0 }, /* 7 */
85  {"Liege", -infinity, -6.385, 30.0, 66.2, 0.0 }, /* 8 */
86  {"Loenhout", 0.0, 4.8, 0.0, 77.0, 2.28 }, /* 9 */
87  {"Mons", -infinity, -6.848, 0.0, 66.2, 0.0 }, /* 10 */
88  {"Namur", -infinity, -2.120, 0.0, 66.2, 0.0 }, /* 11 */
89  {"Petange", -infinity, -1.919, 25.0, 66.2, 0.0 }, /* 12 */
90  {"Peronnes", 0.0, 0.96, 0.0, 66.2, 1.68 }, /* 13 */
91  {"Sinsin", 0.0, 0.0, 0.0, 63.0, 0.0 }, /* 14 */
92  {"Voeren", 20.344, 22.012, 50.0, 66.2, 1.68 }, /* 15 */
93  {"Wanze", 0.0, 0.0, 0.0, 66.2, 0.0 }, /* 16 */
94  {"Warnand", 0.0, 0.0, 0.0, 66.2, 0.0 }, /* 17 */
95  {"Zeebrugge", 8.87, 11.594, 0.0, 77.0, 2.28 }, /* 18 */
96  {"Zomergem", 0.0, 0.0, 0.0, 80.0, 0.0 } /* 19 */
97 };
98 
99 /** data of arcs */
100 static const ArcData arcdata[] =
101 {
102  /* node1 node2 diameter length active */
103  { 18, 6, 890.0, 4.0, FALSE },
104  { 18, 6, 890.0, 4.0, FALSE },
105  { 6, 5, 890.0, 6.0, FALSE },
106  { 6, 5, 890.0, 6.0, FALSE },
107  { 5, 19, 890.0, 26.0, FALSE },
108  { 9, 1, 590.1, 43.0, FALSE },
109  { 1, 7, 590.1, 29.0, FALSE },
110  { 7, 19, 590.1, 19.0, FALSE },
111  { 19, 13, 890.0, 55.0, FALSE },
112  { 15, 3, 890.0, 5.0, TRUE },
113  { 15, 3, 395.0, 5.0, TRUE },
114  { 3, 8, 890.0, 20.0, FALSE },
115  { 3, 8, 395.0, 20.0, FALSE },
116  { 8, 17, 890.0, 25.0, FALSE },
117  { 8, 17, 395.0, 25.0, FALSE },
118  { 17, 11, 890.0, 42.0, FALSE },
119  { 11, 0, 890.0, 40.0, FALSE },
120  { 0, 13, 890.0, 5.0, FALSE },
121  { 13, 10, 890.0, 10.0, FALSE },
122  { 10, 4, 890.0, 25.0, FALSE },
123  { 17, 16, 395.5, 10.5, FALSE },
124  { 16, 14, 315.5, 26.0, TRUE },
125  { 14, 2, 315.5, 98.0, FALSE },
126  { 2, 12, 315.5, 6.0, FALSE }
127 };
128 
129 /** gas temperatur (K) */
130 static const SCIP_Real gastemp = 281.15;
131 
132 /** absolute rugosity (mm) */
133 static const SCIP_Real rugosity = 0.05;
134 
135 /** density of gas relative to air */
136 static const SCIP_Real density = 0.616;
137 
138 /** compressibility factor */
139 static const SCIP_Real compressibility = 0.8;
140 
141 
142 
143 /** sets up problem */
144 static
146  SCIP* scip /**< SCIP data structure */
147  )
148 {
149  SCIP_VAR* flow[narcs];
150  SCIP_VAR* supply[nnodes];
151  SCIP_VAR* pressure[nnodes];
152  SCIP_VAR* pressurediff[narcs];
153 
154  SCIP_CONS* flowbalance[nnodes];
155  SCIP_CONS* pressurediffcons[narcs];
156  SCIP_CONS* pressureloss[narcs];
157 
158  char name[SCIP_MAXSTRLEN];
159  int i;
160  int j;
161 
162  /* create empty problem */
163  SCIP_CALL( SCIPcreateProbBasic(scip, "gastrans") );
164 
165  /* create variables for flows and add to problem */
166  for( i = 0; i < narcs; ++i )
167  {
168  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "flow_%s_%s", nodedata[arcdata[i].node1].name, nodedata[arcdata[i].node2].name);
169  SCIP_CALL( SCIPcreateVarBasic(scip, &flow[i], name, arcdata[i].active ? 0.0 : -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
170 
171  SCIP_CALL( SCIPaddVar(scip, flow[i]) );
172  }
173 
174  /* create variables for pressure difference and add to problem */
175  for( i = 0; i < narcs; ++i )
176  {
177  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "pressurediff_%s_%s", nodedata[arcdata[i].node1].name, nodedata[arcdata[i].node2].name);
178  SCIP_CALL( SCIPcreateVarBasic(scip, &pressurediff[i], name, -SCIPinfinity(scip), SCIPinfinity(scip), 0.0, SCIP_VARTYPE_CONTINUOUS) );
179 
180  SCIP_CALL( SCIPaddVar(scip, pressurediff[i]) );
181  }
182 
183  /* create variables for supply at nodes and add to problem */
184  for( i = 0; i < nnodes; ++i )
185  {
186  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "supply_%s", nodedata[i].name);
187  SCIP_CALL( SCIPcreateVarBasic(scip, &supply[i], name,
188  nodedata[i].supplylower == -infinity ? -SCIPinfinity(scip) : nodedata[i].supplylower,
189  nodedata[i].supplyupper == infinity ? SCIPinfinity(scip) : nodedata[i].supplyupper,
190  nodedata[i].cost, SCIP_VARTYPE_CONTINUOUS) ); /*lint !e777*/
191 
192  SCIP_CALL( SCIPaddVar(scip, supply[i]) );
193  }
194 
195  /* create variables for squared pressure at nodes and add to problem */
196  for( i = 0; i < nnodes; ++i )
197  {
198  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "pressure_%s", nodedata[i].name);
199  SCIP_CALL( SCIPcreateVarBasic(scip, &pressure[i], name,
200  nodedata[i].pressurelower*nodedata[i].pressurelower, nodedata[i].pressureupper*nodedata[i].pressureupper,
201  0.0, SCIP_VARTYPE_CONTINUOUS) );
202 
203  SCIP_CALL( SCIPaddVar(scip, pressure[i]) );
204  }
205 
206  /* create flow balance constraints and add to problem
207  * for each node i: outflows - inflows = supply
208  */
209  for( i = 0; i < nnodes; ++i )
210  {
211  SCIP_Real minusone;
212 
213  minusone = -1.0;
214 
215  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "flowbalance%s", nodedata[i].name);
216  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &flowbalance[i], name, 1, &supply[i], &minusone, 0.0, 0.0) );
217 
218  for( j = 0; j < narcs; ++j )
219  {
220  /* check for outflows */
221  if( arcdata[j].node1 == i )
222  {
223  SCIP_CALL( SCIPaddCoefLinear(scip, flowbalance[i], flow[j], +1.0) );
224  }
225 
226  /* check for inflows */
227  if( arcdata[j].node2 == i )
228  {
229  SCIP_CALL( SCIPaddCoefLinear(scip, flowbalance[i], flow[j], -1.0) );
230  }
231  }
232 
233  SCIP_CALL( SCIPaddCons(scip, flowbalance[i]) );
234  }
235 
236  /* create pressure difference constraints and add to problem
237  * pressurediff[node1 to node2] = pressure[node1] - pressure[node2]
238  */
239  for( i = 0; i < narcs; ++i )
240  {
241  SCIP_VAR* vars[3];
242  SCIP_Real coefs[3] = { 1.0, -1.0, 1.0 };
243 
244  assert(pressurediff[i] != NULL);
245  assert(arcdata[i].node1 >= 0 && arcdata[i].node1 < nnodes);
246  assert(arcdata[i].node2 >= 0 && arcdata[i].node2 < nnodes);
247  assert(pressure[arcdata[i].node1] != NULL);
248  assert(pressure[arcdata[i].node2] != NULL);
249 
250  vars[0] = pressurediff[i];
251  vars[1] = pressure[arcdata[i].node1];
252  vars[2] = pressure[arcdata[i].node2];
253 
254  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "pressurediff_%s_%s", nodedata[arcdata[i].node1].name, nodedata[arcdata[i].node2].name);
255  SCIP_CALL( SCIPcreateConsBasicLinear(scip, &pressurediffcons[i], name, 3, vars, coefs, 0.0, 0.0) );
256 
257  SCIP_CALL( SCIPaddCons(scip, pressurediffcons[i]) );
258  }
259 
260  /* create pressure loss constraints and add to problem
261  * pressure loss on active arcs: flow(i) * flow(i) + coef * pressurediff(i) <= 0.0
262  * on regular pipes: flow(i) * abs(flow(i)) - coef * pressurediff(i) == 0.0,
263  * where coef = 96.074830e-15*power(diameter(i)^5/lambda/compressibility/temperatur/length(i)/density
264  * and lambda = (2*log10(3.7*diameter(i)/rugosity))^(-2);
265  */
266  for( i = 0; i < narcs; ++i )
267  {
268  SCIP_EXPR* exprflow;
269  SCIP_EXPR* exprs[2];
270  SCIP_Real coefs[2];
271  SCIP_EXPR* exprsum;
272  SCIP_Real coef;
273 
274  coef = 96.074830e-15 * pow(arcdata[i].diameter, 5.0) * pow(2.0*log10(3.7*arcdata[i].diameter / rugosity), 2.0)
275  / compressibility / gastemp / arcdata[i].length / density;
276 
277  /* we can always use the signpower-expression, because flow(i) is positive for active arcs */
278  SCIP_CALL( SCIPcreateExprVar(scip, &exprflow, flow[i], NULL, NULL) );
279  SCIP_CALL( SCIPcreateExprSignpower(scip, &exprs[0], exprflow, 2.0, NULL, NULL) );
280  SCIP_CALL( SCIPreleaseExpr(scip, &exprflow) );
281 
282  SCIP_CALL( SCIPcreateExprVar(scip, &exprs[1], pressurediff[i], NULL, NULL) );
283 
284  coefs[0] = 1.0;
285  coefs[1] = arcdata[i].active ? coef : -coef;
286  SCIP_CALL( SCIPcreateExprSum(scip, &exprsum, 2, exprs, coefs, 0.0, NULL, NULL) );
287  SCIP_CALL( SCIPreleaseExpr(scip, &exprs[1]) );
288  SCIP_CALL( SCIPreleaseExpr(scip, &exprs[0]) );
289 
290  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "pressureloss_%s_%s", nodedata[arcdata[i].node1].name, nodedata[arcdata[i].node2].name);
291  SCIP_CALL( SCIPcreateConsBasicNonlinear(scip, &pressureloss[i], name, exprsum, arcdata[i].active ? -SCIPinfinity(scip) : 0.0, 0.0) );
292  SCIP_CALL( SCIPreleaseExpr(scip, &exprsum) );
293 
294  SCIP_CALL( SCIPaddCons(scip, pressureloss[i]) );
295  }
296 
297  /* release variables and constraints
298  * the problem has them captured, and we do not require them anymore
299  */
300  for( i = 0; i < nnodes; ++i )
301  {
302  SCIP_CALL( SCIPreleaseVar(scip, &supply[i]) );
303  SCIP_CALL( SCIPreleaseVar(scip, &pressure[i]) );
304 
305  SCIP_CALL( SCIPreleaseCons(scip, &flowbalance[i]) );
306  }
307 
308  for( i = 0; i < narcs; ++i )
309  {
310  SCIP_CALL( SCIPreleaseVar(scip, &flow[i]) );
311  SCIP_CALL( SCIPreleaseVar(scip, &pressurediff[i]) );
312 
313  SCIP_CALL( SCIPreleaseCons(scip, &pressurediffcons[i]) );
314  SCIP_CALL( SCIPreleaseCons(scip, &pressureloss[i]) );
315  }
316 
317  return SCIP_OKAY;
318 }
319 
320 /** runs gas transportation example */
321 static
323 {
324  SCIP* scip;
325 
326  SCIP_CALL( SCIPcreate(&scip) );
328 
329  SCIPinfoMessage(scip, NULL, "\n");
330  SCIPinfoMessage(scip, NULL, "*****************************************\n");
331  SCIPinfoMessage(scip, NULL, "* Running Gas Transportation Model *\n");
332  SCIPinfoMessage(scip, NULL, "*****************************************\n");
333  SCIPinfoMessage(scip, NULL, "\n");
334 
335  SCIP_CALL( setupProblem(scip) );
336 
337  SCIPinfoMessage(scip, NULL, "Original problem:\n");
338  SCIP_CALL( SCIPprintOrigProblem(scip, NULL, "cip", FALSE) );
339 
340  SCIPinfoMessage(scip, NULL, "\n");
341  SCIP_CALL( SCIPpresolve(scip) );
342 
343  /* SCIPinfoMessage(scip, NULL, "Presolved problem:\n");
344  SCIP_CALL( SCIPprintTransProblem(scip, NULL, "cip", FALSE) );
345  */
346 
347  SCIPinfoMessage(scip, NULL, "\nSolving...\n");
348  SCIP_CALL( SCIPsolve(scip) );
349 
350  /*
351  if( SCIPgetNSols(scip) > 0 )
352  {
353  SCIPinfoMessage(scip, NULL, "\nSolution:\n");
354  SCIP_CALL( SCIPprintSol(scip, SCIPgetBestSol(scip), NULL, FALSE) );
355  }
356  */
357 
358  SCIP_CALL( SCIPfree(&scip) );
359 
360  return SCIP_OKAY;
361 }
362 
363 
364 /** main method starting SCIP */
365 int main(
366  int argc, /**< number of arguments from the shell */
367  char** argv /**< array of shell arguments */
368  )
369 { /*lint --e{715}*/
370  SCIP_RETCODE retcode;
371 
372  retcode = runGastrans();
373 
374  /* evaluate return code of the SCIP process */
375  if( retcode != SCIP_OKAY )
376  {
377  /* write error back trace */
378  SCIPprintError(retcode);
379  return -1;
380  }
381 
382  return 0;
383 }
#define narcs
Definition: gastrans.c:68
SCIP_Real cost
Definition: gastrans.c:52
SCIP_Real supplyupper
Definition: gastrans.c:49
SCIP_RETCODE SCIPcreateConsBasicLinear(SCIP *scip, SCIP_CONS **cons, const char *name, int nvars, SCIP_VAR **vars, SCIP_Real *vals, SCIP_Real lhs, SCIP_Real rhs)
#define infinity
Definition: gastrans.c:71
SCIP_Real pressureupper
Definition: gastrans.c:51
#define SCIP_MAXSTRLEN
Definition: def.h:293
SCIP_RETCODE SCIPcreateExprSignpower(SCIP *scip, SCIP_EXPR **expr, SCIP_EXPR *child, SCIP_Real exponent, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: expr_pow.c:3190
SCIP_Real supplylower
Definition: gastrans.c:48
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip_var.c:1245
static SCIP_RETCODE runGastrans(void)
Definition: gastrans.c:322
#define FALSE
Definition: def.h:87
SCIP_Real pressurelower
Definition: gastrans.c:50
SCIP_Real length
Definition: gastrans.c:60
SCIP_Real SCIPinfinity(SCIP *scip)
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10755
#define TRUE
Definition: def.h:86
static const SCIP_Real gastemp
Definition: gastrans.c:130
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:54
SCIP_RETCODE SCIPcreateVarBasic(SCIP *scip, SCIP_VAR **var, const char *name, SCIP_Real lb, SCIP_Real ub, SCIP_Real obj, SCIP_VARTYPE vartype)
Definition: scip_var.c:185
static GRAPHNODE ** active
struct NodeData NodeData
SCIP_RETCODE SCIPcreateExprVar(SCIP *scip, SCIP_EXPR **expr, SCIP_VAR *var, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: expr_var.c:377
SCIP_RETCODE SCIPcreate(SCIP **scip)
Definition: scip_general.c:283
SCIP_RETCODE SCIPcreateExprSum(SCIP *scip, SCIP_EXPR **expr, int nchildren, SCIP_EXPR **children, SCIP_Real *coefficients, SCIP_Real constant, SCIP_DECL_EXPR_OWNERCREATE((*ownercreate)), void *ownercreatedata)
Definition: expr_sum.c:1070
SCIP_RETCODE SCIPaddCoefLinear(SCIP *scip, SCIP_CONS *cons, SCIP_VAR *var, SCIP_Real val)
void SCIPinfoMessage(SCIP *scip, FILE *file, const char *formatstr,...)
Definition: scip_message.c:199
SCIP_RETCODE SCIPcreateProbBasic(SCIP *scip, const char *name)
Definition: scip_prob.c:170
static const NodeData nodedata[]
Definition: gastrans.c:74
SCIP_RETCODE SCIPsolve(SCIP *scip)
Definition: scip_solve.c:2613
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip_prob.c:2768
SCIP_Bool active
Definition: gastrans.c:61
static const SCIP_Real rugosity
Definition: gastrans.c:133
SCIP_RETCODE SCIPcreateConsBasicNonlinear(SCIP *scip, SCIP_CONS **cons, const char *name, SCIP_EXPR *expr, SCIP_Real lhs, SCIP_Real rhs)
SCIP_RETCODE SCIPpresolve(SCIP *scip)
Definition: scip_solve.c:2443
SCIP_RETCODE SCIPprintOrigProblem(SCIP *scip, FILE *file, const char *extension, SCIP_Bool genericnames)
#define NULL
Definition: lpi_spx1.cpp:155
#define SCIP_CALL(x)
Definition: def.h:384
SCIP_Real diameter
Definition: gastrans.c:59
int node2
Definition: gastrans.c:58
#define SCIP_Bool
Definition: def.h:84
SCIP_RETCODE SCIPincludeDefaultPlugins(SCIP *scip)
SCIP_RETCODE SCIPreleaseExpr(SCIP *scip, SCIP_EXPR **expr)
Definition: scip_expr.c:1407
struct ArcData ArcData
static const ArcData arcdata[]
Definition: gastrans.c:100
static SCIP_RETCODE setupProblem(SCIP *scip)
Definition: gastrans.c:145
const char * name
Definition: gastrans.c:47
static const SCIP_Real density
Definition: gastrans.c:136
int node1
Definition: gastrans.c:57
int main(int argc, char **argv)
Definition: gastrans.c:365
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
Definition: scip_prob.c:1666
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
Definition: scip_cons.c:1110
#define SCIP_Real
Definition: def.h:177
#define nnodes
Definition: gastrans.c:65
void SCIPprintError(SCIP_RETCODE retcode)
Definition: scip_general.c:211
default SCIP plugins
SCIP callable library.
SCIP_RETCODE SCIPfree(SCIP **scip)
Definition: scip_general.c:315
static const SCIP_Real compressibility
Definition: gastrans.c:139