1 // $Id$
2 #include "flopc.hpp"
3 using namespace flopc;
4 using namespace std;
5 #include <OsiCbcSolverInterface.hpp>
6
7 // M A G I C Power Scheduling Problem
8 /* A number of power stations are committed to meet demand for a particular
9 day. three types of generators having different operating characteristics
10 are available. Generating units can be shut down or operate between
11 minimum and maximum output levels. Units can be started up or closed down
12 in every demand block.
13
14 Reference: R E Day and H P Williams, MAGIC: The design and use of an
15 interactive modeling language for mathematical programming,
16 Department Business Studies, University of Edinburgh, Edinburgh,
17 UK, February 1982.
18
19 H P Williams, Model Building in Mathematical Programming, Wiley,
20 New York, 1978.
21 */
22
main()23 int main() {
24 enum {t12pm_6am, t6am_9am, t9am_3pm, t3pm_6pm, t6pm_12pm, numT};
25 enum {type_1, type_2, type_3, numG};
26
27 MP_set T(numT); T.cyclic();
28 MP_set G(numG);
29
30 MP_data dem(T); // demand (1000mw)
31 dem(t12pm_6am) = 15;
32 dem(t6am_9am) = 30;
33 dem(t9am_3pm) = 25;
34 dem(t3pm_6pm) = 40;
35 dem(t6pm_12pm) = 27;
36
37 MP_data dur(T); // duration (hours)
38 dur(t12pm_6am) = 6;
39 dur(t6am_9am) = 3;
40 dur(t9am_3pm) = 6;
41 dur(t3pm_6pm) = 3;
42 dur(t6pm_12pm) = 6;
43
44 MP_data
45 min_pow(G),
46 max_pow(G),
47 cost_min(G),
48 cost_inc(G),
49 start(G),
50 number(G);
51
52 min_pow(type_1) = 0.85; max_pow(type_1) = 2.0;
53 min_pow(type_2) = 1.25; max_pow(type_2) = 1.75;
54 min_pow(type_3) = 1.5; max_pow(type_3) = 4.0;
55
56 cost_min(type_1) = 1000; cost_inc(type_1) = 2.0;
57 cost_min(type_2) = 2600; cost_inc(type_2) = 1.3;
58 cost_min(type_3) = 3000; cost_inc(type_3) = 3.0;
59
60 start(type_1) = 2000; number(type_1) = 12;
61 start(type_2) = 1000; number(type_2) = 10;
62 start(type_3) = 500; number(type_3) = 5;
63
64 MP_data
65 peak, // peak power (1000 mw)
66 ener(T), // energy demand in load block (1000mwh)
67 tener, // total energy demanded (1000mwh)
68 lf; // load factor
69
70 MP_index i,j;
71
72 peak() = maximum(T(i), dem(i));
73
74 ener(i) = dur(i)*dem(i);
75
76 tener() = sum(T(i), ener(i));
77
78 lf() = tener()/(peak() * 24);
79
80 peak.display();
81 tener.display();
82 lf.display();
83 ener.display("ener");
84
85 MP_variable
86 x(G,T), // generator output (1000mw)
87 n(G,T), // number of generators in use
88 s(G,T); // number of generators started up
89
90 n.integer();
91 n.upperLimit(j,i) = number(j);
92 n.upperLimit.display("n upper");
93
94 MP_constraint
95 pow(T), // demand for power (1000mw)
96 res(T), // spinning reserve requirements (1000mw)
97 st(G,T), // start_up definition
98 minu(G,T), // minimum generation level (1000mw)
99 maxu(G,T); // maximum generation level (1000mw)
100
101 pow(i) = sum(G(j), x(j,i)) >= dem(i);
102
103 res(i) = sum(G(j), max_pow(j)*n(j,i)) >= 1.15*dem(i);
104
105 st(j,i) = s(j,i) >= n(j,i) - n(j,i-1);
106
107 minu(j,i) = x(j,i) >= min_pow(j)*n(j,i);
108
109 maxu(j,i) = x(j,i) <= max_pow(j)*n(j,i);
110
111 MP_expression cost; // total operating cost (l)
112 cost = sum(G(j)*T(i),
113 dur(i)* cost_min(j)*n(j,i) + start(j)*s(j,i) +
114 1000*dur(i)*cost_inc(j)*(x(j,i)-min_pow(j)*n(j,i)) );
115
116
117 MP_model::getDefaultModel().setSolver(new OsiCbcSolverInterface);
118 MP_model::getDefaultModel().verbose();
119
120 minimize(cost);
121
122 assert(MP_model::getDefaultModel()->getNumRows()==55);
123 assert(MP_model::getDefaultModel()->getNumCols()==45);
124 assert(MP_model::getDefaultModel()->getNumElements()==135);
125 assert(MP_model::getDefaultModel()->getObjValue()>=988539 && MP_model::getDefaultModel()->getObjValue()<=988541);
126
127
128 x.display("x");
129 n.display("n");
130 s.display("s");
131 const int TSize=T.size();
132 double **rep = new double *[TSize];
133 for(int cnt=0;cnt<TSize;cnt++)
134 {
135 rep[cnt]=new double[4];
136 }
137 for (unsigned i=0; i<T.size(); i++) {
138 rep[i][0] = dem(i);
139
140 rep[i][1] = 0;
141 for (unsigned j=0; j<G.size(); j++) {
142 rep[i][1] += max_pow(j) * n.level(j,i);
143 }
144
145 rep[i][2] = 0;
146 for (unsigned j=0; j<G.size(); j++) {
147 rep[i][2] += s.level(j,i);
148 }
149
150 rep[i][3] = -pow.price(i)/dur(i)/1000;
151 }
152
153 for (unsigned i=0; i<T.size(); i++) {
154 cout<<i<<" "<<rep[i][0]<<" "<<rep[i][1]<<" "<<rep[i][2]<<" "<<rep[i][3]<<endl;
155 }
156
157 cout<<"Test magic passed."<<endl;
158 }
159