CONOPT
Loading...
Searching...
No Matches
pinadd.c
Go to the documentation of this file.
1
6
7#include <stdlib.h>
8#include <stdio.h>
9#include <math.h>
10#include <string.h>
11#include "coiheader.h"
12
13#include "comdecl.h"
14int T = 16; /* Number of time periods */
15#define Tmin 16
16#define Tmax 20
17#define Vpp 7 /* Variables per period */
18#define Epp 6 /* Equations per period */
19double Xkeep[Tmax*Vpp]; /* Space for storing solutions */
20int Xstat[Tmax*Vpp]; /* Space for storing status info */
21int Estat[Tmax*Epp]; /* Space for storing status info */
22
23
28int COI_CALLCONV Pin_ReadMatrix( double LOWER[], double CURR[], double UPPER[], int VSTA[], int TYPE[],
29 double RHS[], int ESTA[], int COLSTA[], int ROWNO[], double VALUE[],
30 int NLFLAG[], int NUMVAR, int NUMCON, int NUMNZ, void* USRMEM )
31{
32 int it, is, i, icol, iz, iold;
33/*
34 Define the information for the columns.
35
36 We should not supply status information, VSTA.
37
38 We order the variables as follows:
39 0: td, 1: cs, 2: s, 3: d, 4: r, 5: p, and 6: rev. All variables for
40 period 1 appears first followed by all variables for period 2 etc.
41
42 td, cs, s, and d have lower bounds of 0, r and p have lower
43 bounds of 1, and rev has no lower bound.
44 All have infinite upper bounds (default).
45 The initial value of td is 18, s is 7, cs is 7*t, d is td-s,
46 p is 14, and r is r(t-1)-d. No initial value for rev. */
47
48 for ( it=1; it<=T; it++ ) {
49 is = Vpp*(it-1);
50 LOWER[is ] = 0;
51 LOWER[is+1] = 0;
52 LOWER[is+2] = 0;
53 LOWER[is+3] = 0;
54 LOWER[is+4] = 1;
55 LOWER[is+5] = 1;
56 CURR[is ] = 18;
57 CURR[is+1] = 7*it;
58 CURR[is+2] = 7;
59 CURR[is+3] = CURR[is] - CURR[is+2];
60 if ( it > 1 )
61 CURR[is+4] = CURR[is+4-Vpp] - CURR[is+3];
62 else
63 CURR[is+4] = 500 - CURR[is+3];
64 CURR[is+5] = 14;
65 }
66 if ( T > Tmin ) {
67
68/* This is a restart: Use the initial values from last solve for
69 the variables in the first periods and extrapolate the last
70 period using a linear extrapolation between the last two */
71
72 iold = Vpp*(T-1);
73 for ( i = 0; i < iold; i++ )
74 CURR[i] = Xkeep[i];
75 for ( i = 0; i < Vpp; i++ )
76 CURR[iold+i] = CURR[iold+i-Vpp] + (CURR[iold+i-Vpp]-CURR[iold+i-2*Vpp]);
77
78/* The variables from the last solve are given the old status and
79 the variables in the new period are given those in the last
80 pariod. Similarly with the Equation status: */
81
82 for ( i = 0; i < iold; i++ )
83 VSTA[i] = Xstat[i];
84 for ( i = 0; i < Vpp; i++ )
85 VSTA[iold+i] = Xstat[iold+i-Vpp];
86 iold = Epp*(T-1)+1;
87 for ( i = 0; i < iold; i++ )
88 ESTA[i] = Estat[i];
89 for ( i = 0; i < Epp; i++ )
90 ESTA[iold+i] = Estat[iold+i-Epp];
91 }
92/*
93 Define the information for the rows
94
95 We order the constraints as follows: The objective is first,
96 followed by tddef, sdef, csdef, ddef, rdef, and revdef for
97 the first period, the same for the second period, etc.
98
99 The objective is a nonbinding constraint: */
100
101 TYPE[0] = 3;
102
103/* All others are equalities: */
104
105 for ( i = 1; i < NUMCON; i++ )
106 TYPE[i] = 0;
107/*
108 Right hand sides: In all periods except the first, only tddef
109 has a nonzero right hand side of 1+2.3*1.015**(t-1).
110 In the initial period there are contributions from lagged
111 variables in the constraints that have lagged variables. */
112
113 for ( it = 1; it<=T; it++ ) {
114 is = 1 + 6*(it-1);
115 RHS[is] = 1.0+2.3*pow(1.015,(it-1.));
116 }
117
118/* tddef: + 0.87*td(0) */
119
120 RHS[1] = RHS[1]+ 0.87*18.;
121
122/* sdef: +0.75*s(0) */
123
124 RHS[2] = 0.75*6.5;
125
126/* csdef: +1*cs(0) */
127
128 RHS[3] = 0;
129
130/* rdef: +1*r(0) */
131
132 RHS[5] = 500;
133
134/*
135 Define the structure and content of the Jacobian:
136 To help define the Jacobian pattern and values it can be useful to
137 make a picture of the Jacobian. We describe the variables for one
138 period and the constraints they are part of:
139
140 td cs s d r p rev
141 0 1 2 3 4 5 6
142 0: Obj (1+r)**(1-t)
143 Period t:
144 0: tddef 1.0 0.13
145 1: sdef NL 1.0 NL
146 2: csdef 1.0 -1.0
147 3: ddef -1.0 1.0 1.0
148 4: rdef 1.0 1.0
149 5: revdef NL NL NL 1.0
150 Period t+1:
151 0: tddef -0.87
152 1: sdef -0.75
153 2: csdef -1.0
154 3: ddef
155 4: rdef -1.0
156 5: revdef
157
158 The Jacobian has to be sorted column-wise so we will just define
159 the elements column by column according to the table above: */
160
161 iz = 0; /* Running index for Jacobian in C convention 0 to nz-1 */
162 icol = 0; /* Running column index in C convention 0 to N-1 */
163 for ( it = 1; it <= T; it++ )
164 {
165
166/* is points to the position before the first equation for the period */
167
168 is = 1 + 6*(it-1); /* 1 correspond to the objective in row 0 */
169
170/* Column td: */
171 COLSTA[icol] = iz;
172 icol = icol + 1;
173 ROWNO[iz] = is;
174 VALUE[iz] = +1.0;
175 NLFLAG[iz] = 0;
176 iz = iz + 1;
177 ROWNO[iz] = is+3;
178 VALUE[iz] = -1.0;
179 NLFLAG[iz] = 0;
180 iz = iz + 1;
181 if ( it < T ) {
182 ROWNO[iz] = is+Epp;
183 VALUE[iz] = -0.87;
184 NLFLAG[iz] = 0;
185 iz = iz + 1;
186 }
187
188/* Column cs */
189 COLSTA[icol] = iz;
190 icol = icol + 1;
191 ROWNO[iz] = is+1;
192 NLFLAG[iz] = 1;
193 iz = iz + 1;
194 ROWNO[iz] = is+2;
195 VALUE[iz] = +1.0;
196 NLFLAG[iz] = 0;
197 iz = iz + 1;
198 if ( it < T ) {
199 ROWNO[iz] = is+2+Epp;
200 VALUE[iz] = -1.0;
201 NLFLAG[iz] = 0;
202 iz = iz + 1;
203 }
204
205/* Column s */
206 COLSTA[icol] = iz;
207 icol = icol + 1;
208 ROWNO[iz] = is+1;
209 VALUE[iz] = +1.0;
210 NLFLAG[iz] = 0;
211 iz = iz + 1;
212 ROWNO[iz] = is+2;
213 VALUE[iz] = -1.0;
214 NLFLAG[iz] = 0;
215 iz = iz + 1;
216 ROWNO[iz] = is+3;
217 VALUE[iz] = +1.0;
218 NLFLAG[iz] = 0;
219 iz = iz + 1;
220 if ( it < T ) {
221 ROWNO[iz] = is+1+Epp;
222 VALUE[iz] = -0.75;
223 NLFLAG[iz] = 0;
224 iz = iz + 1;
225 }
226
227/* Column d: */
228 COLSTA[icol] = iz;
229 icol = icol + 1;
230 ROWNO[iz] = is+3;
231 VALUE[iz] = +1.0;
232 NLFLAG[iz] = 0;
233 iz = iz + 1;
234 ROWNO[iz] = is+4;
235 VALUE[iz] = +1.0;
236 NLFLAG[iz] = 0;
237 iz = iz + 1;
238 ROWNO[iz] = is+5;
239 NLFLAG[iz] = 1;
240 iz = iz + 1;
241
242/* Column r: */
243 COLSTA[icol] = iz;
244 icol = icol + 1;
245 ROWNO[iz] = is+4;
246 VALUE[iz] = +1.0;
247 NLFLAG[iz] = 0;
248 iz = iz + 1;
249 ROWNO[iz] = is+5;
250 NLFLAG[iz] = 1;
251 iz = iz + 1;
252 if ( it < T ) {
253 ROWNO[iz] = is+4+Epp;
254 VALUE[iz] = -1.0;
255 NLFLAG[iz] = 0;
256 iz = iz + 1;
257 }
258
259/* Column p: */
260 COLSTA[icol] = iz;
261 icol = icol + 1;
262 ROWNO[iz] = is;
263 VALUE[iz] = +0.13;
264 NLFLAG[iz] = 0;
265 iz = iz + 1;
266 ROWNO[iz] = is+1;
267 NLFLAG[iz] = 1;
268 iz = iz + 1;
269 ROWNO[iz] = is+5;
270 NLFLAG[iz] = 1;
271 iz = iz + 1;
272
273/* Column rev: */
274 COLSTA[icol] = iz;
275 icol = icol + 1;
276 ROWNO[iz] = +0; /* Objective -- no is term */
277 VALUE[iz] = pow(1.05,(1.0-it));
278 NLFLAG[iz] = 0;
279 iz = iz + 1;
280 ROWNO[iz] = is+5;
281 VALUE[iz] = 1.0;
282 NLFLAG[iz] = 0;
283 iz = iz + 1;
284 }
285 COLSTA[icol] = iz; /* First elements past last column */
286
287 return 0;
288}
289/*
290 =====================================================================
291 Compute nonlinear terms and non-constant Jacobian elements */
292
293
298int COI_CALLCONV Pin_FDEval( const double X[], double* G, double JAC[], int ROWNO, const int JACNUM[], int MODE,
299 int IGNERR, int* ERRCNT, int NUMVAR, int NUMJAC, int THREAD, void* USRMEM )
300{
301
302 int it, is;
303 double h1, h2;
304
305/* Compute the number of the period, it */
306
307 it = (ROWNO+Epp-1) / Epp;
308 is = Vpp*(it-1);
309 if ( 1+(it-1)*Epp+1 == ROWNO ) {
310
311/* sdef equation. Nonlinear term = -(1.1+0.1*p)*1.02**(-cs/7) */
312
313 h1 = (1.1+0.1*X[is+5]);
314 h2 = pow(1.02,-X[is+1]/7.0);
315 if ( 1 == MODE || 3 == MODE )
316 *G = -h1*h2;
317 if ( 2 == MODE || 3 == MODE ) {
318 JAC[is+1] = h1*h2*log(1.02)/7.0;
319 JAC[is+5] = -h2*0.1;
320 }
321 }
322 else if ( 1+(it-1)*Epp+5 == ROWNO ) {
323
324/* revdef equation. Nonlinear term = -d*(p-250/r) */
325
326 if ( 1 == MODE || 3 == MODE )
327 *G = -X[is+3]*(X[is+5]-250./X[is+4]);
328 if ( 2 == MODE || 3 == MODE ) {
329 JAC[is+3] = -(X[is+5]-250./X[is+4]);
330 JAC[is+4] = -X[is+3]*250./pow(X[is+4],2);
331 JAC[is+5] = -X[is+3];
332 }
333 }
334 else {
335
336/* Error - this equation is not nonlinear */
337
338 printf("\nError. Pin_FDEval called with ROWNO = %d. MODE = %d\n\n", ROWNO, MODE);
339 return 1;
340 };
341 return 0;
342}
343
344int COI_CALLCONV Pin_Solution( const double XVAL[], const double XMAR[], const int XBAS[], const int XSTA[],
345 const double YVAL[], const double YMAR[], const int YBAS[], const int YSTA[],
346 int NUMVAR, int NUMCON, void* USRMEM )
347{
348/* Specialized Solution routine */
349 int i;
350 char *status[4] = {"Lower","Upper","Basic","Super"};
351 if ( T < Tmax ) {
352 printf("Saving primal solution and status information for T = %d\n",T);
353 for ( i = 0; i < NUMVAR; i++ ) {
354 Xkeep[i] = XVAL[i];
355 Xstat[i] = XBAS[i];
356 }
357 for ( i = 0; i < NUMCON; i++ )
358 Estat[i] = YBAS[i];
359 return 0;
360 }
361 fprintf(fd,"\n Variable Solution value Reduced cost Status\n\n");
362 for ( i=0; i<NUMVAR; i++ )
363 fprintf(fd,"%6d%18f%18f%10s\n", i, XVAL[i], XMAR[i], status[XBAS[i]] );
364 fprintf(fd,"\n Constrnt Activity level Marginal cost Status\n\n");
365 for ( i=0; i<NUMCON; i++ )
366 fprintf(fd,"%6d%18f%18f%10s\n", i, YVAL[i], YMAR[i], status[YBAS[i]] );
367
368 return 0;
369}
370
371#include "std.c"
372
375int main(int argc, char** argv)
376{
377 c_log( "Starting to execute", START );
378
379 T = Tmin;
380
381/* Tell CONOPT about the sizes in the model */
382/* Number of variables (excl. slacks): 7 per period */
383
384 COI_Error += COIDEF_NumVar ( CntVect, 7*T );
385
386/* Number of equations: 1 objective + 6 per period */
387
388 COI_Error += COIDEF_NumCon ( CntVect, 1+6*T );
389
390/* Number of nonzeros in the Jacobian. See the counting in */
391/* ReadMatrix above. For each period there is 1 in the objective, 16 */
392/* for unlagged variables and 4 for lagged variables. */
393
394 COI_Error += COIDEF_NumNz ( CntVect, 17*T+4*(T-1) );
395
396/* Number of nonlinear nonzeros. 5 unlagged for each period. */
397
398 COI_Error += COIDEF_NumNlNz ( CntVect, 5*T );
399
400/* Objective: Maximize Constraint no 0 */
401
402 COI_Error += COIDEF_OptDir ( CntVect, 1 );
403 COI_Error += COIDEF_ObjCon ( CntVect, 0 );
404
405/* Turn debugging of FDEval on or off */
406
407 COI_Error += COIDEF_DebugFV ( CntVect, 0 );
408
409/* Register the options file */
410
411 COI_Error += COIDEF_Optfile ( CntVect, "pindyck.opt");
412 /*
413 Register the necessary callback routines with CONOPT
414 */
415 COI_Error += COIDEF_Message ( CntVect, &Std_Message ); /* Register the callback Message */
416 COI_Error += COIDEF_ErrMsg ( CntVect, &Std_ErrMsg ); /* Register the callback ErrMsg */
417 COI_Error += COIDEF_Status ( CntVect, &Std_Status ); /* Register the callback Status */
418 COI_Error += COIDEF_Solution ( CntVect, &Pin_Solution ); /* Register the callback Solution */
419 COI_Error += COIDEF_ReadMatrix( CntVect, &Pin_ReadMatrix); /* Register the callback ReadMatrix */
420 COI_Error += COIDEF_FDEval ( CntVect, &Pin_FDEval); /* Register the callback FDEval */
421
422#if defined(LICENSE_INT_1) && defined(LICENSE_INT_2) && defined(LICENSE_INT_3) && defined(LICENSE_TEXT)
423 COI_Error += COIDEF_License ( CntVect, LICENSE_INT_1, LICENSE_INT_2, LICENSE_INT_3, LICENSE_TEXT);
424#endif
425
426 if ( COI_Error ) {
427 printf("Skipping COI_Solve due to setup errors. COI_Error = %d\n",COI_Error);
428 c_log( "Skipping Solve due to setup errors", COI_Error);
429 }
430 COI_Error = COI_Solve ( CntVect ); /* Optimize */
431 printf("After solving. COI_Error =%d\n",COI_Error);
432 if ( COI_Error ) {
433 c_log( "Errors encountered during first solve", COI_Error);
434 }
435 else if ( mstat != 2 || sstat != 1 ) {
436 c_log( "Incorrect Model or Solver Status during first solve", -1);
437 }
438
439/* Now increase the time horizon gradually up to 20.
440 All sizes that depend on T must be registered again. */
441
442 for ( T = Tmin+1; T <= Tmax; T++ ) {
443 COI_Error += COIDEF_NumVar ( CntVect, 7*T );
444 COI_Error += COIDEF_NumCon ( CntVect, 1+6*T );
445 COI_Error += COIDEF_NumNz ( CntVect, 17*T+4*(T-1) );
446 COI_Error += COIDEF_NumNlNz ( CntVect, 5*T );
447
448/* This time we have status information copied from a previous solve,
449 i.e. IniStat = 2 */
450
451 COI_Error += COIDEF_IniStat ( CntVect, 2 );
452 if ( COI_Error ) {
453 printf("\n**** Fatal Error while revising CONOPT Sizes.\n");
454 c_log( "Errors encountered while revising model sizes", COI_Error);
455 }
456 COI_Error = COI_Solve( CntVect );
457 if ( COI_Error ) {
458 printf("\n**** Fatal Error while Solving for T=%d\n",T);
459 c_log( "Errors encountered during later solve", COI_Error);
460 }
461 else if ( mstat != 2 || sstat != 1 ) {
462 printf("\n**** Incorrect Model or Solver Status returned for T=%d\n",T);
463 c_log( "Incorrect Model or Solver Status during later solve", -1);
464 }
465 }
466 printf("\nEnd of Pinadd Model. Return code=%d\n",COI_Error);
467 c_log( "Successful Solve", OK );
468}
coiheader.h
C language header file for direct linking against the COI library generated by apiwrapper for GAMS Ve...
CntVect
coiHandle_t CntVect
Definition comdecl.h:14
START
#define START
Definition comdecl.h:13
fd
FILE * fd
Definition comdecl.h:3
mstat
int mstat
Definition comdecl.h:7
OK
#define OK
Definition comdecl.h:12
sstat
int sstat
Definition comdecl.h:8
COI_Error
int COI_Error
Definition comdecl.h:15
COIDEF_ReadMatrix
int COI_CALLCONV COIDEF_ReadMatrix(coiHandle_t cntvect, COI_READMATRIX_t coi_readmatrix)
define callback routine for providing the matrix data to CONOPT.
COIDEF_Message
int COI_CALLCONV COIDEF_Message(coiHandle_t cntvect, COI_MESSAGE_t coi_message)
define callback routine for handling messages returned during the solution process.
COIDEF_Solution
int COI_CALLCONV COIDEF_Solution(coiHandle_t cntvect, COI_SOLUTION_t coi_solution)
define callback routine for returning the final solution values.
COIDEF_Status
int COI_CALLCONV COIDEF_Status(coiHandle_t cntvect, COI_STATUS_t coi_status)
define callback routine for returning the completion status.
COIDEF_ErrMsg
int COI_CALLCONV COIDEF_ErrMsg(coiHandle_t cntvect, COI_ERRMSG_t coi_errmsg)
define callback routine for returning error messages for row, column or Jacobian elements.
COIDEF_FDEval
int COI_CALLCONV COIDEF_FDEval(coiHandle_t cntvect, COI_FDEVAL_t coi_fdeval)
define callback routine for performing function and derivative evaluations.
COIDEF_Optfile
int COI_CALLCONV COIDEF_Optfile(coiHandle_t cntvect, const char *optfile)
define callback routine for defining an options file.
COIDEF_DebugFV
int COI_CALLCONV COIDEF_DebugFV(coiHandle_t cntvect, int debugfv)
turn Debugging of FDEval on and off.
COIDEF_License
int COI_CALLCONV COIDEF_License(coiHandle_t cntvect, int licint1, int licint2, int licint3, const char *licstring)
define the License Information.
COIDEF_IniStat
int COI_CALLCONV COIDEF_IniStat(coiHandle_t cntvect, int inistat)
handling of the initial status values.
COIDEF_ObjCon
int COI_CALLCONV COIDEF_ObjCon(coiHandle_t cntvect, int objcon)
defines the Objective Constraint.
COIDEF_NumVar
int COI_CALLCONV COIDEF_NumVar(coiHandle_t cntvect, int numvar)
defines the number of variables in the model.
COIDEF_NumNz
int COI_CALLCONV COIDEF_NumNz(coiHandle_t cntvect, int numnz)
defines the number of nonzero elements in the Jacobian.
COIDEF_NumCon
int COI_CALLCONV COIDEF_NumCon(coiHandle_t cntvect, int numcon)
defines the number of constraints in the model.
COIDEF_OptDir
int COI_CALLCONV COIDEF_OptDir(coiHandle_t cntvect, int optdir)
defines the Optimization Direction.
COIDEF_NumNlNz
int COI_CALLCONV COIDEF_NumNlNz(coiHandle_t cntvect, int numnlnz)
defines the Number of Nonlinear Nonzeros.
COI_Solve
int COI_CALLCONV COI_Solve(coiHandle_t cntvect)
method for starting the solving process of CONOPT.
Pin_ReadMatrix
int COI_CALLCONV Pin_ReadMatrix(double LOWER[], double CURR[], double UPPER[], int VSTA[], int TYPE[], double RHS[], int ESTA[], int COLSTA[], int ROWNO[], double VALUE[], int NLFLAG[], int NUMVAR, int NUMCON, int NUMNZ, void *USRMEM)
Define information about the model.
Definition pinadd.c:28
main
int main(int argc, char **argv)
Main program. A simple setup and call of CONOPT.
Definition pinadd.c:375
Xstat
int Xstat[Tmax *Vpp]
Definition pinadd.c:20
Pin_Solution
int COI_CALLCONV Pin_Solution(const double XVAL[], const double XMAR[], const int XBAS[], const int XSTA[], const double YVAL[], const double YMAR[], const int YBAS[], const int YSTA[], int NUMVAR, int NUMCON, void *USRMEM)
Definition pinadd.c:344
Xkeep
double Xkeep[Tmax *Vpp]
Definition pinadd.c:19
Pin_FDEval
int COI_CALLCONV Pin_FDEval(const double X[], double *G, double JAC[], int ROWNO, const int JACNUM[], int MODE, int IGNERR, int *ERRCNT, int NUMVAR, int NUMJAC, int THREAD, void *USRMEM)
Compute nonlinear terms and non-constant Jacobian elements.
Definition pinadd.c:298
T
int T
Definition pinadd.c:14
Tmin
#define Tmin
Definition pinadd.c:15
Estat
int Estat[Tmax *Epp]
Definition pinadd.c:21
Epp
#define Epp
Definition pinadd.c:18
Tmax
#define Tmax
Definition pinadd.c:16
Vpp
#define Vpp
Definition pinadd.c:17
Std_Status
int COI_CALLCONV Std_Status(int MODSTA, int SOLSTA, int ITER, double OBJVAL, void *USRMEM)
Definition std.c:45
Std_Message
int COI_CALLCONV Std_Message(int SMSG, int DMSG, int NMSG, char *MSGV[], void *USRMEM)
Definition std.c:8
c_log
void c_log(char *msgt, int code)
Definition std.c:202
Std_ErrMsg
int COI_CALLCONV Std_ErrMsg(int ROWNO, int COLNO, int POSNO, const char *MSG, void *USRMEM)
Definition std.c:26