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-2017 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 email to scip@zib.de. */
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 abspower constraints in SCIP 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  * http://www.gams.com/modlib/libhtml/gastrans.htm
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[2]
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_Real coef;
269 
270  (void) SCIPsnprintf(name, SCIP_MAXSTRLEN, "pressureloss_%s_%s", nodedata[arcdata[i].node1].name, nodedata[arcdata[i].node2].name);
271 
272  coef = 96.074830e-15 * pow(arcdata[i].diameter, 5.0) * pow(2.0*log10(3.7*arcdata[i].diameter / rugosity), 2.0)
273  / compressibility / gastemp / arcdata[i].length / density;
274 
275  if( arcdata[i].active )
276  {
277  /* we can also use an abspower constraint here, because flow(i) is positive for active arcs */
278  SCIP_CALL( SCIPcreateConsBasicAbspower(scip, &pressureloss[i], name, flow[i], pressurediff[i], 2.0, 0.0, coef, -SCIPinfinity(scip), 0.0) );
279  }
280  else
281  {
282  SCIP_CALL( SCIPcreateConsBasicAbspower(scip, &pressureloss[i], name, flow[i], pressurediff[i], 2.0, 0.0, -coef, 0.0, 0.0) );
283  }
284 
285  SCIP_CALL( SCIPaddCons(scip, pressureloss[i]) );
286  }
287 
288  /* release variables and constraints
289  * the problem has them captured, and we do not require them anymore
290  */
291  for( i = 0; i < nnodes; ++i )
292  {
293  SCIP_CALL( SCIPreleaseVar(scip, &supply[i]) );
294  SCIP_CALL( SCIPreleaseVar(scip, &pressure[i]) );
295 
296  SCIP_CALL( SCIPreleaseCons(scip, &flowbalance[i]) );
297  }
298 
299  for( i = 0; i < narcs; ++i )
300  {
301  SCIP_CALL( SCIPreleaseVar(scip, &flow[i]) );
302  SCIP_CALL( SCIPreleaseVar(scip, &pressurediff[i]) );
303 
304  SCIP_CALL( SCIPreleaseCons(scip, &pressurediffcons[i]) );
305  SCIP_CALL( SCIPreleaseCons(scip, &pressureloss[i]) );
306  }
307 
308  return SCIP_OKAY;
309 }
310 
311 /* runs gas transportation example */
312 static
314 {
315  SCIP* scip;
316 
317  SCIP_CALL( SCIPcreate(&scip) );
319 
320  SCIPinfoMessage(scip, NULL, "\n");
321  SCIPinfoMessage(scip, NULL, "*****************************************\n");
322  SCIPinfoMessage(scip, NULL, "* Running Gas Transportation Model *\n");
323  SCIPinfoMessage(scip, NULL, "*****************************************\n");
324  SCIPinfoMessage(scip, NULL, "\n");
325 
326  SCIP_CALL( setupProblem(scip) );
327 
328  SCIPinfoMessage(scip, NULL, "Original problem:\n");
329  SCIP_CALL( SCIPprintOrigProblem(scip, NULL, "cip", FALSE) );
330 
331  SCIPinfoMessage(scip, NULL, "\n");
332  SCIP_CALL( SCIPpresolve(scip) );
333 
334  /* SCIPinfoMessage(scip, NULL, "Presolved problem:\n");
335  SCIP_CALL( SCIPprintTransProblem(scip, NULL, "cip", FALSE) );
336  */
337 
338  SCIPinfoMessage(scip, NULL, "\nSolving...\n");
339  SCIP_CALL( SCIPsolve(scip) );
340 
341  SCIP_CALL( SCIPfreeTransform(scip) );
342 
343  /*
344  if( SCIPgetNSols(scip) > 0 )
345  {
346  SCIPinfoMessage(scip, NULL, "\nSolution:\n");
347  SCIP_CALL( SCIPprintSol(scip, SCIPgetBestSol(scip), NULL, FALSE) );
348  }
349  */
350 
351  SCIP_CALL( SCIPfree(&scip) );
352 
353  return SCIP_OKAY;
354 }
355 
356 
357 /** main method starting SCIP */
358 int main(
359  int argc, /**< number of arguments from the shell */
360  char** argv /**< array of shell arguments */
361  ) /*lint --e{715}*/
362 {
363  SCIP_RETCODE retcode;
364 
365  retcode = runGastrans();
366 
367  /* evaluate return code of the SCIP process */
368  if( retcode != SCIP_OKAY )
369  {
370  /* write error back trace */
371  SCIPprintError(retcode);
372  return -1;
373  }
374 
375  return 0;
376 }
#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:259
SCIPInterval pow(const SCIPInterval &x, const SCIPInterval &y)
SCIP_Real supplylower
Definition: gastrans.c:48
SCIP_RETCODE SCIPreleaseVar(SCIP *scip, SCIP_VAR **var)
Definition: scip.c:18760
static SCIP_RETCODE runGastrans(void)
Definition: gastrans.c:313
#define FALSE
Definition: def.h:64
SCIP_Real pressurelower
Definition: gastrans.c:50
SCIP_Real length
Definition: gastrans.c:60
SCIP_Real SCIPinfinity(SCIP *scip)
Definition: scip.c:47022
int SCIPsnprintf(char *t, int len, const char *s,...)
Definition: misc.c:10011
#define TRUE
Definition: def.h:63
static const SCIP_Real gastemp
Definition: gastrans.c:130
enum SCIP_Retcode SCIP_RETCODE
Definition: type_retcode.h:53
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.c:17693
static GRAPHNODE ** active
struct NodeData NodeData
SCIP_RETCODE SCIPcreate(SCIP **scip)
Definition: scip.c:748
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.c:1343
SCIP_RETCODE SCIPcreateProbBasic(SCIP *scip, const char *name)
Definition: scip.c:9997
SCIP_RETCODE SCIPcreateConsBasicAbspower(SCIP *scip, SCIP_CONS **cons, const char *name, SCIP_VAR *x, SCIP_VAR *z, SCIP_Real exponent, SCIP_Real xoffset, SCIP_Real zcoef, SCIP_Real lhs, SCIP_Real rhs)
static const NodeData nodedata[]
Definition: gastrans.c:74
SCIP_RETCODE SCIPsolve(SCIP *scip)
Definition: scip.c:16109
SCIP_RETCODE SCIPaddCons(SCIP *scip, SCIP_CONS *cons)
Definition: scip.c:12585
SCIP_Bool active
Definition: gastrans.c:61
static const SCIP_Real rugosity
Definition: gastrans.c:133
SCIP_RETCODE SCIPpresolve(SCIP *scip)
Definition: scip.c:15948
SCIP_RETCODE SCIPprintOrigProblem(SCIP *scip, FILE *file, const char *extension, SCIP_Bool genericnames)
Definition: scip.c:44212
#define SCIP_CALL(x)
Definition: def.h:350
SCIP_Real diameter
Definition: gastrans.c:59
int node2
Definition: gastrans.c:58
SCIP_RETCODE SCIPfreeTransform(SCIP *scip)
Definition: scip.c:17246
#define SCIP_Bool
Definition: def.h:61
SCIP_RETCODE SCIPincludeDefaultPlugins(SCIP *scip)
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:358
SCIP_RETCODE SCIPaddVar(SCIP *scip, SCIP_VAR *var)
Definition: scip.c:11486
SCIP_RETCODE SCIPreleaseCons(SCIP *scip, SCIP_CONS **cons)
Definition: scip.c:27755
#define SCIP_Real
Definition: def.h:149
#define nnodes
Definition: gastrans.c:65
void SCIPprintError(SCIP_RETCODE retcode)
Definition: scip.c:676
default SCIP plugins
SCIP callable library.
SCIP_RETCODE SCIPfree(SCIP **scip)
Definition: scip.c:780
static const SCIP_Real compressibility
Definition: gastrans.c:139