trunk/models/pipeline.a4c

(*
	ASCEND modelling environment
	Copyright (C) 1998, 2007 Carnegie Mellon University

	This program is free software; you can redistribute it and/or modify
	it under the terms of the GNU General Public License as published by
	the Free Software Foundation; either version 2, or (at your option)
	any later version.

	This program is distributed in the hope that it will be useful,
	but WITHOUT ANY WARRANTY; without even the implied warranty of
	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
	GNU General Public License for more details.
	
	You should have received a copy of the GNU General Public License
	along with this program; if not, write to the Free Software
	Foundation, Inc., 59 Temple Place - Suite 330,
	Boston, MA 02111-1307, USA.
*)
REQUIRE "atoms.a4l";
PROVIDE "pipeline.a4c";
(*
	The following the conditional model is discussed in the PhD thesis of
	Vicente Rico-Ramirez,
	https://pse.cheme.cmu.edu/ascend/ftp/pdfThesis/victhesis.pdf

	The problem consists of a  network of 38 pipes, 5 of them including a check
	valve (so that the flow in those pipes is possible only in one direction).
	Each pipe is represented by a set of algebraic disjunctive equations. It
	represents a problem which we can represent as a conditional model and 
	solve with the ASCEND conditional solver CMSlv.

	The problem was originally described in the paper by Bullard and Biegler,
	
	Bullard and Biegler, Iterated Linear Programming Strategies for Nonsmooth
	Simulation: Continuous and Mixed Integer Approaches. Comp. Chem. Eng.
	Vol. 16, 949-961, 1992

	Note that you will need CONOPT installed on your system in order to solve
	this problem.

	Model by Vicente Rico-Ramirez, April 10, 1998
*)

(*---------------------------
	Base model of a pipe
*)
MODEL arc;
	k		IS_A k_constant;
	Q		IS_A volume_rate;
	H		IS_A distance;
METHODS
END arc;

(*---------------------------
	Pipe without check valve
*)
MODEL arc_w_no_valve REFINES arc;
	bol		IS_A boolean_var;

	(* Boundary *)
	CONDITIONAL
		cond: Q >= 0.0 {gpm};
	END CONDITIONAL;

	bol == SATISFIED(cond,1e-08 {gpm});

	(* Variant Equations *)
	eq1: H = k * sqr(Q);
	eq2: H = -k * sqr(Q);

	(* Disjunctive Statement *)
    WHEN (bol)
	CASE TRUE:
			USE eq1;
	CASE FALSE:
			USE eq2;
	END WHEN;

METHODS
    METHOD default_self;
    END default_self;

    METHOD specify;
    	FIX k;
    END specify;

    METHOD bound_self;
		Q.lower_bound := -1e50{gpm};
		H.lower_bound := -1e50{ft};
    END bound_self;

    METHOD values;
		H := 0 {ft};
		bol := SATISFIED(cond,1e-08 {gpm});
    END values;

END arc_w_no_valve;

(*---------------------------
	Pipe with check valve
*)
MODEL arc_w_valve REFINES arc;

	bol		IS_A boolean_var;

	(* Boundary *)
	CONDITIONAL
		cond: H >= 0.0 {ft};
	END CONDITIONAL ;

	bol == SATISFIED(cond,1e-08{ft});

	(* Variant Equations *)
	eq1: sqr(Q) = H/k;
	eq2: Q = 0.0 {gpm};

	(* Disjunctive Statement *)
    WHEN (bol)
	CASE TRUE:
		USE eq1;
	CASE FALSE:
		USE eq2;
	END WHEN;

METHODS
    METHOD default_self;
    END default_self;

    METHOD specify;
    	FIX k;
    END specify;

    METHOD bound_self;
		H.lower_bound := -1e50{ft};
    END bound_self;

    METHOD values;
		H := 1.0 {ft};
		bol := SATISFIED(cond,1e-08{ft});
    END values;
END arc_w_valve;

(*---------------------------
	Pipeline Network 
*)
MODEL pipeline;

	nodes		IS_A set OF integer_constant;
	arcs		IS_A set OF integer_constant;
	Pipe[arcs]	IS_A arc;
	P[nodes]	IS_A distance;
	W[nodes]	IS_A volume_rate;
	Q[arcs]		IS_A volume_rate;
	H[arcs]		IS_A distance;
	k		IS_A k_constant;

	(* Sets definition *)
	nodes :== [1..22];
	arcs  :== [1..38];

	(* Now wire-up the network *)
	Pipe[15],Pipe[19],Pipe[23],Pipe[32],
	    Pipe[34] IS_REFINED_TO arc_w_valve;

	FOR i IN [1..14] CREATE
	    Pipe[i] IS_REFINED_TO arc_w_no_valve;
	END FOR;
	FOR i IN [16..18] CREATE
	    Pipe[i] IS_REFINED_TO arc_w_no_valve;
	END FOR;
	FOR i IN [20..22] CREATE
	    Pipe[i] IS_REFINED_TO arc_w_no_valve;
	END FOR ;
	FOR i IN [24..31] CREATE
	    Pipe[i] IS_REFINED_TO arc_w_no_valve;
	END FOR;
	Pipe[33] IS_REFINED_TO arc_w_no_valve;
	FOR i IN [35..38] CREATE
	    Pipe[i] IS_REFINED_TO arc_w_no_valve;
	END FOR;

	FOR i IN arcs CREATE
	    Q[i],Pipe[i].Q ARE_THE_SAME;
	    H[i],Pipe[i].H ARE_THE_SAME;
	    k,Pipe[i].k ARE_THE_SAME;
	END FOR;

	(* Invariant Equations *)

	(* Pressure Drops at each arc *)
	P[1]  - P[6]  = H[1];
	P[2]  - P[1]  = H[2];
	P[2]  - P[3]  = H[3];
	P[3]  - P[11] = H[4];
	P[13] - P[11] = H[5];
	P[14] - P[12] = H[6];
	P[12] - P[16] = H[7];
	P[16] - P[17] = H[8];
	P[19] - P[20] = H[9];
	P[21] - P[19] = H[10];
	P[7]  - P[6]  = H[11];
	P[10] - P[9]  = H[12];
	P[7]  - P[9]  = H[13];
	P[21] - P[18] = H[14];
	P[14] - P[10] = H[15];
	P[22] - P[10] = H[16];
	P[21] - P[22] = H[17];
	P[1]  - P[5]  = H[18];
	P[5]  - P[4]  = H[19];
	P[4]  - P[2]  = H[20];
	P[3]  - P[4]  = H[21];
	P[11] - P[4]  = H[22];
	P[11] - P[12] = H[23];
	P[12] - P[13] = H[24];
	P[13] - P[14] = H[25];
	P[12] - P[15] = H[26];
	P[16] - P[15] = H[27];
	P[18] - P[15] = H[28];
	P[17] - P[18] = H[29];
	P[19] - P[15] = H[30];
	P[15] - P[20] = H[31];
	P[20] - P[21] = H[32];
	P[9]  - P[22] = H[33];
	P[9]  - P[8]  = H[34];
	P[5]  - P[8]  = H[35];
	P[5]  - P[10] = H[36];
	P[8]  - P[7]  = H[37];
	P[6]  - P[5]  = H[38];

	(* Flow balance around each node  *)

	W[1] = Q[2] - Q[1] - Q[18];
	W[2] = Q[20] - Q[2] - Q[3];
	W[3] = Q[3] - Q[4] - Q[21];
	W[4] = Q[19] + Q[21] + Q[22] - Q[20];
	W[5] = Q[18] + Q[38] - Q[19] - Q[35] - Q[36];
	W[6] = Q[1] + Q[11] - Q[38];
	W[7] = Q[37] - Q[11] - Q[13];
	W[8] = Q[34] + Q[35] - Q[37];
	W[9] = Q[12] + Q[13] - Q[33] - Q[34];
	W[10] = Q[15] + Q[16] + Q[36] - Q[12];
	W[11] = Q[4] + Q[5] - Q[22] - Q[23];
	W[12] = Q[6] + Q[23] - Q[7] - Q[24] - Q[26];
	W[13] = Q[24] - Q[5] - Q[25];
	W[14] = Q[25] - Q[6] - Q[15];
	W[15] = Q[26] + Q[27] + Q[28] + Q[30] - Q[31];
	W[16] = Q[7] - Q[8] - Q[27];
	W[17] = Q[8] - Q[29];
	W[18] = Q[14] + Q[29] - Q[28];
	W[19] = Q[10] - Q[9] - Q[30];
	W[20] = Q[9] + Q[31] - Q[32];
	W[21] = Q[32] - Q[10] - Q[14] - Q[17];
	W[22] = Q[17] + Q[33] - Q[16];

METHODS
    METHOD default_self;
    END default_self;

    METHOD specify;
		FIX k;
		FOR i IN nodes DO
		    FIX W[i];
		END FOR;
		FREE W[18];
		FIX P[17];
    END specify;

	METHOD bound_self;
		FOR i IN nodes DO
		    W[i].lower_bound := -1e50 {gpm};
		END FOR;
		FOR i IN nodes DO
		    P[i].lower_bound := 0 {ft};
		END FOR ;
	END bound_self;

	METHOD bound_all;
		FOR i IN arcs DO
		    RUN Pipe[i].bound_self;
		END FOR;
	    RUN bound_self;
	END bound_all;

    METHOD values;
		RUN bound_all;
		FOR i IN arcs DO
			RUN Pipe[i].values;
		END FOR;

		(* fixed *)
		k := 0.36412197 {s^2/ft^5};

		FOR i IN nodes DO
			W[i] := 0 {gpm};
		END FOR;
		W[1] := -897.6 {gpm};
		W[7] := -1570.9 {gpm};
		W[11] := 897.6 {gpm};
		W[17] := 1795.3 {gpm};
		W[20] := 448.8 {gpm};
		W[22] := -673.2 {gpm};

		(* Initial guess *)

		(* Pressure head *)
		FOR i IN nodes DO
			P[i] := 0 {ft};
		END FOR ;
		P[17] := 0 {ft};

		(* volume rates *)
		Q[1] :=  -48.5 {gpm};
		Q[2] :=  -640.7 {gpm};
		Q[3] :=  393.2 {gpm};
		Q[4] :=  538.9 {gpm};
		Q[5] :=  -32.3 {gpm};
		Q[6] :=  490.2 {gpm};
		Q[7] :=  649.9 {gpm};
		Q[8] :=  904.7 {gpm};
		Q[9] :=  112.2 {gpm};
		Q[10] :=  232.5 {gpm};
		Q[11] :=  402.3 {gpm};
		Q[12] :=  -420.4 {gpm};
		Q[13] :=  687.3 {gpm};
		Q[14] :=  621.7 {gpm};
		Q[15] :=  -719.2 {gpm};
		Q[16] :=  52.7 {gpm};
		Q[17] :=  1199 {gpm};
		Q[18] :=  305.4 {gpm};
		Q[19] :=  582.8 {gpm};
		Q[20] :=  -247.5 {gpm};
		Q[21] :=  -145.7 {gpm};
		Q[22] :=  -684.6 {gpm};
		Q[23] :=  293.6 {gpm};
		Q[24] :=  -261.3 {gpm};
		Q[25] :=  -229 {gpm};
		Q[26] :=  -395.1 {gpm};
		Q[27] :=  -254.8 {gpm};
		Q[28] :=  -268.9 {gpm};
		Q[29] :=  -890.6 {gpm};
		Q[30] :=  120.3 {gpm};
		Q[31] :=  -8.1{gpm};
		Q[32] :=  -344.7 {gpm};
		Q[33] :=  473.1 {gpm};
		Q[34] :=  -206.2 {gpm};
		Q[35] :=  -275 {gpm};
		Q[36] :=  351.5 {gpm};
		Q[37] :=  -481.2 {gpm};
		Q[38] :=  353.9 {gpm};
	END values;

	METHOD on_load;
		RUN default_self;
		RUN reset;
		RUN values;
	END on_load;

	METHOD self_test;
		(* no tests yet *)
	END self_test;

END pipeline;
(* :ex: set ts=4: *)
1	(*
2	ASCEND modelling environment
3	Copyright (C) 1998, 2007 Carnegie Mellon University
4
5	This program is free software; you can redistribute it and/or modify
6	it under the terms of the GNU General Public License as published by
7	the Free Software Foundation; either version 2, or (at your option)
8	any later version.
9
10	This program is distributed in the hope that it will be useful,
11	but WITHOUT ANY WARRANTY; without even the implied warranty of
12	MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13	GNU General Public License for more details.
14
15	You should have received a copy of the GNU General Public License
16	along with this program; if not, write to the Free Software
17	Foundation, Inc., 59 Temple Place - Suite 330,
18	Boston, MA 02111-1307, USA.
19	*)
20	REQUIRE "atoms.a4l";
21	PROVIDE "pipeline.a4c";
22	(*
23	The following the conditional model is discussed in the PhD thesis of
24	Vicente Rico-Ramirez,
25	https://pse.cheme.cmu.edu/ascend/ftp/pdfThesis/victhesis.pdf
26
27	The problem consists of a network of 38 pipes, 5 of them including a check
28	valve (so that the flow in those pipes is possible only in one direction).
29	Each pipe is represented by a set of algebraic disjunctive equations. It
30	represents a problem which we can represent as a conditional model and
31	solve with the ASCEND conditional solver CMSlv.
32
33	The problem was originally described in the paper by Bullard and Biegler,
34
35	Bullard and Biegler, Iterated Linear Programming Strategies for Nonsmooth
36	Simulation: Continuous and Mixed Integer Approaches. Comp. Chem. Eng.
37	Vol. 16, 949-961, 1992
38
39	Note that you will need CONOPT installed on your system in order to solve
40	this problem.
41
42	Model by Vicente Rico-Ramirez, April 10, 1998
43	*)
44
45	(*---------------------------
46	Base model of a pipe
47	*)
48	MODEL arc;
49	k IS_A k_constant;
50	Q IS_A volume_rate;
51	H IS_A distance;
52	METHODS
53	END arc;
54
55	(*---------------------------
56	Pipe without check valve
57	*)
58	MODEL arc_w_no_valve REFINES arc;
59	bol IS_A boolean_var;
60
61	(* Boundary *)
62	CONDITIONAL
63	cond: Q >= 0.0 {gpm};
64	END CONDITIONAL;
65
66	bol == SATISFIED(cond,1e-08 {gpm});
67
68	(* Variant Equations *)
69	eq1: H = k * sqr(Q);
70	eq2: H = -k * sqr(Q);
71
72	(* Disjunctive Statement *)
73	WHEN (bol)
74	CASE TRUE:
75	USE eq1;
76	CASE FALSE:
77	USE eq2;
78	END WHEN;
79
80	METHODS
81	METHOD default_self;
82	END default_self;
83
84	METHOD specify;
85	FIX k;
86	END specify;
87
88	METHOD bound_self;
89	Q.lower_bound := -1e50{gpm};
90	H.lower_bound := -1e50{ft};
91	END bound_self;
92
93	METHOD values;
94	H := 0 {ft};
95	bol := SATISFIED(cond,1e-08 {gpm});
96	END values;
97
98	END arc_w_no_valve;
99
100	(*---------------------------
101	Pipe with check valve
102	*)
103	MODEL arc_w_valve REFINES arc;
104
105	bol IS_A boolean_var;
106
107	(* Boundary *)
108	CONDITIONAL
109	cond: H >= 0.0 {ft};
110	END CONDITIONAL ;
111
112	bol == SATISFIED(cond,1e-08{ft});
113
114	(* Variant Equations *)
115	eq1: sqr(Q) = H/k;
116	eq2: Q = 0.0 {gpm};
117
118	(* Disjunctive Statement *)
119	WHEN (bol)
120	CASE TRUE:
121	USE eq1;
122	CASE FALSE:
123	USE eq2;
124	END WHEN;
125
126	METHODS
127	METHOD default_self;
128	END default_self;
129
130	METHOD specify;
131	FIX k;
132	END specify;
133
134	METHOD bound_self;
135	H.lower_bound := -1e50{ft};
136	END bound_self;
137
138	METHOD values;
139	H := 1.0 {ft};
140	bol := SATISFIED(cond,1e-08{ft});
141	END values;
142	END arc_w_valve;
143
144	(*---------------------------
145	Pipeline Network
146	*)
147	MODEL pipeline;
148
149	nodes IS_A set OF integer_constant;
150	arcs IS_A set OF integer_constant;
151	Pipe[arcs] IS_A arc;
152	P[nodes] IS_A distance;
153	W[nodes] IS_A volume_rate;
154	Q[arcs] IS_A volume_rate;
155	H[arcs] IS_A distance;
156	k IS_A k_constant;
157
158	(* Sets definition *)
159	nodes :== [1..22];
160	arcs :== [1..38];
161
162	(* Now wire-up the network *)
163	Pipe[15],Pipe[19],Pipe[23],Pipe[32],
164	Pipe[34] IS_REFINED_TO arc_w_valve;
165
166	FOR i IN [1..14] CREATE
167	Pipe[i] IS_REFINED_TO arc_w_no_valve;
168	END FOR;
169	FOR i IN [16..18] CREATE
170	Pipe[i] IS_REFINED_TO arc_w_no_valve;
171	END FOR;
172	FOR i IN [20..22] CREATE
173	Pipe[i] IS_REFINED_TO arc_w_no_valve;
174	END FOR ;
175	FOR i IN [24..31] CREATE
176	Pipe[i] IS_REFINED_TO arc_w_no_valve;
177	END FOR;
178	Pipe[33] IS_REFINED_TO arc_w_no_valve;
179	FOR i IN [35..38] CREATE
180	Pipe[i] IS_REFINED_TO arc_w_no_valve;
181	END FOR;
182
183	FOR i IN arcs CREATE
184	Q[i],Pipe[i].Q ARE_THE_SAME;
185	H[i],Pipe[i].H ARE_THE_SAME;
186	k,Pipe[i].k ARE_THE_SAME;
187	END FOR;
188
189	(* Invariant Equations *)
190
191	(* Pressure Drops at each arc *)
192	P[1] - P[6] = H[1];
193	P[2] - P[1] = H[2];
194	P[2] - P[3] = H[3];
195	P[3] - P[11] = H[4];
196	P[13] - P[11] = H[5];
197	P[14] - P[12] = H[6];
198	P[12] - P[16] = H[7];
199	P[16] - P[17] = H[8];
200	P[19] - P[20] = H[9];
201	P[21] - P[19] = H[10];
202	P[7] - P[6] = H[11];
203	P[10] - P[9] = H[12];
204	P[7] - P[9] = H[13];
205	P[21] - P[18] = H[14];
206	P[14] - P[10] = H[15];
207	P[22] - P[10] = H[16];
208	P[21] - P[22] = H[17];
209	P[1] - P[5] = H[18];
210	P[5] - P[4] = H[19];
211	P[4] - P[2] = H[20];
212	P[3] - P[4] = H[21];
213	P[11] - P[4] = H[22];
214	P[11] - P[12] = H[23];
215	P[12] - P[13] = H[24];
216	P[13] - P[14] = H[25];
217	P[12] - P[15] = H[26];
218	P[16] - P[15] = H[27];
219	P[18] - P[15] = H[28];
220	P[17] - P[18] = H[29];
221	P[19] - P[15] = H[30];
222	P[15] - P[20] = H[31];
223	P[20] - P[21] = H[32];
224	P[9] - P[22] = H[33];
225	P[9] - P[8] = H[34];
226	P[5] - P[8] = H[35];
227	P[5] - P[10] = H[36];
228	P[8] - P[7] = H[37];
229	P[6] - P[5] = H[38];
230
231	(* Flow balance around each node *)
232
233	W[1] = Q[2] - Q[1] - Q[18];
234	W[2] = Q[20] - Q[2] - Q[3];
235	W[3] = Q[3] - Q[4] - Q[21];
236	W[4] = Q[19] + Q[21] + Q[22] - Q[20];
237	W[5] = Q[18] + Q[38] - Q[19] - Q[35] - Q[36];
238	W[6] = Q[1] + Q[11] - Q[38];
239	W[7] = Q[37] - Q[11] - Q[13];
240	W[8] = Q[34] + Q[35] - Q[37];
241	W[9] = Q[12] + Q[13] - Q[33] - Q[34];
242	W[10] = Q[15] + Q[16] + Q[36] - Q[12];
243	W[11] = Q[4] + Q[5] - Q[22] - Q[23];
244	W[12] = Q[6] + Q[23] - Q[7] - Q[24] - Q[26];
245	W[13] = Q[24] - Q[5] - Q[25];
246	W[14] = Q[25] - Q[6] - Q[15];
247	W[15] = Q[26] + Q[27] + Q[28] + Q[30] - Q[31];
248	W[16] = Q[7] - Q[8] - Q[27];
249	W[17] = Q[8] - Q[29];
250	W[18] = Q[14] + Q[29] - Q[28];
251	W[19] = Q[10] - Q[9] - Q[30];
252	W[20] = Q[9] + Q[31] - Q[32];
253	W[21] = Q[32] - Q[10] - Q[14] - Q[17];
254	W[22] = Q[17] + Q[33] - Q[16];
255
256	METHODS
257	METHOD default_self;
258	END default_self;
259
260	METHOD specify;
261	FIX k;
262	FOR i IN nodes DO
263	FIX W[i];
264	END FOR;
265	FREE W[18];
266	FIX P[17];
267	END specify;
268
269	METHOD bound_self;
270	FOR i IN nodes DO
271	W[i].lower_bound := -1e50 {gpm};
272	END FOR;
273	FOR i IN nodes DO
274	P[i].lower_bound := 0 {ft};
275	END FOR ;
276	END bound_self;
277
278	METHOD bound_all;
279	FOR i IN arcs DO
280	RUN Pipe[i].bound_self;
281	END FOR;
282	RUN bound_self;
283	END bound_all;
284
285	METHOD values;
286	RUN bound_all;
287	FOR i IN arcs DO
288	RUN Pipe[i].values;
289	END FOR;
290
291	(* fixed *)
292	k := 0.36412197 {s^2/ft^5};
293
294	FOR i IN nodes DO
295	W[i] := 0 {gpm};
296	END FOR;
297	W[1] := -897.6 {gpm};
298	W[7] := -1570.9 {gpm};
299	W[11] := 897.6 {gpm};
300	W[17] := 1795.3 {gpm};
301	W[20] := 448.8 {gpm};
302	W[22] := -673.2 {gpm};
303
304	(* Initial guess *)
305
306	(* Pressure head *)
307	FOR i IN nodes DO
308	P[i] := 0 {ft};
309	END FOR ;
310	P[17] := 0 {ft};
311
312	(* volume rates *)
313	Q[1] := -48.5 {gpm};
314	Q[2] := -640.7 {gpm};
315	Q[3] := 393.2 {gpm};
316	Q[4] := 538.9 {gpm};
317	Q[5] := -32.3 {gpm};
318	Q[6] := 490.2 {gpm};
319	Q[7] := 649.9 {gpm};
320	Q[8] := 904.7 {gpm};
321	Q[9] := 112.2 {gpm};
322	Q[10] := 232.5 {gpm};
323	Q[11] := 402.3 {gpm};
324	Q[12] := -420.4 {gpm};
325	Q[13] := 687.3 {gpm};
326	Q[14] := 621.7 {gpm};
327	Q[15] := -719.2 {gpm};
328	Q[16] := 52.7 {gpm};
329	Q[17] := 1199 {gpm};
330	Q[18] := 305.4 {gpm};
331	Q[19] := 582.8 {gpm};
332	Q[20] := -247.5 {gpm};
333	Q[21] := -145.7 {gpm};
334	Q[22] := -684.6 {gpm};
335	Q[23] := 293.6 {gpm};
336	Q[24] := -261.3 {gpm};
337	Q[25] := -229 {gpm};
338	Q[26] := -395.1 {gpm};
339	Q[27] := -254.8 {gpm};
340	Q[28] := -268.9 {gpm};
341	Q[29] := -890.6 {gpm};
342	Q[30] := 120.3 {gpm};
343	Q[31] := -8.1{gpm};
344	Q[32] := -344.7 {gpm};
345	Q[33] := 473.1 {gpm};
346	Q[34] := -206.2 {gpm};
347	Q[35] := -275 {gpm};
348	Q[36] := 351.5 {gpm};
349	Q[37] := -481.2 {gpm};
350	Q[38] := 353.9 {gpm};
351	END values;
352
353	METHOD on_load;
354	RUN default_self;
355	RUN reset;
356	RUN values;
357	END on_load;
358
359	METHOD self_test;
360	(* no tests yet *)
361	END self_test;
362
363	END pipeline;
364	(* :ex: set ts=4: *)