1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
|
/* For copyright information, see olden_v1.0/COPYRIGHT */
/* main.c
*
* By: Martin C. Carlisle
* 6/16/94
*
* main routine for the Power Pricing problem
* adapted from code by: Steve Lumetta, Sherry Li, and Ismail Khalil
* University of California at Berkeley
*
* Cleaned up the CM stuff. George Necula
*/
#include "power.h"
/* Domain of thetaR->P map is 0.65 to 1.00 [index*0.01+0.65] */
double map_P[36] =
{8752.218091048, 8446.106670416, 8107.990680283,
7776.191574285, 7455.920518777, 7146.602181352,
6847.709026813, 6558.734204024, 6279.213382291,
6008.702199986, 5746.786181029, 5493.078256495,
5247.206333097, 5008.828069358, 4777.615815166,
4553.258735900, 4335.470002316, 4123.971545694,
3918.501939675, 3718.817618538, 3524.683625800,
3335.876573044, 3152.188635673, 2973.421417103,
2799.382330486, 2629.892542617, 2464.782829705,
2303.889031418, 2147.054385395, 1994.132771399,
1844.985347313, 1699.475053321, 1557.474019598,
1418.860479043, 1283.520126656, 1151.338004216};
#define MIN_THETA_R 0.65
#define PER_INDEX_R 0.01
#define MAX_THETA_R 0.995
/* Domain of thetaI->Q map is 0.130 to 0.200 [index*0.002+0.130] */
double map_Q[36] =
{1768.846590190, 1706.229490046, 1637.253873079,
1569.637451623, 1504.419525242, 1441.477913810,
1380.700660446, 1321.980440476, 1265.218982201,
1210.322424636, 1157.203306183, 1105.780028163,
1055.974296746, 1007.714103979, 960.930643875,
915.558722782, 871.538200178, 828.810882006,
787.322098340, 747.020941334, 707.858376214,
669.787829741, 632.765987756, 596.751545633,
561.704466609, 527.587580585, 494.365739051,
462.004890691, 430.472546686, 399.738429196,
369.773787595, 340.550287137, 312.041496095,
284.222260660, 257.068973074, 230.557938283};
#define MIN_THETA_I 0.13
#define PER_INDEX_I 0.002
#define MAX_THETA_I 0.199
#ifdef PLAIN
double wallclock;
#endif
main(int argc,char *argv[])
{
Root r;
int i,finished=0;
double d_theta_R,d_theta_I;
chatting("Past initialization\n");
timer_start(0);
/* initial pass */
r = build_tree();
chatting("Built tree\n");
Compute_Tree(r);
r->last.P = r->D.P;
r->last.Q = r->D.Q;
r->last_theta_R = r->theta_R;
r->last_theta_I = r->theta_I;
r->theta_R = 0.7;
r->theta_I = 0.14;
while (!finished) {
Compute_Tree(r);
chatting("TR=%13.9f, TI=%13.9f, P0=%13.9f, Q0=%13.9f\n",
r->theta_R,r->theta_I,r->D.P,r->D.Q);
if (fabs(r->D.P/10000.0 - r->theta_R) < ROOT_EPSILON &&
fabs(r->D.Q/10000.0 - r->theta_I) < ROOT_EPSILON) {
finished = 1;
} else {
i = (int)((r->theta_R - MIN_THETA_R) / PER_INDEX_R);
if (i<0) i=0;
if (i>35) i=35;
d_theta_R = -(r->theta_R - r->D.P/10000.0) /
(1 - (map_P[i+1] - map_P[i]) / (PER_INDEX_R * 10000.0));
i = (int)((r->theta_I - MIN_THETA_I) / PER_INDEX_I);
if (i<0) i=0;
if (i>35) i=35;
d_theta_I = -(r->theta_I - r->D.Q/10000.0) /
(1 - (map_Q[i+1] - map_Q[i]) / (PER_INDEX_I * 10000.0));
chatting("D TR-%13.9f, TI=%13.9f\n", d_theta_R,d_theta_I);
r->last.P = r->D.P;
r->last.Q = r->D.Q;
r->last_theta_R = r->theta_R;
r->last_theta_I = r->theta_I;
r->theta_R = r->theta_R + d_theta_R;
r->theta_I = r->theta_I + d_theta_I;
}
} /* while */
timer_stop(0);
chatting("Elapsed time %f\n", timer_elapsed(0));
#ifdef FUTURES
__ShutDown(0);
#endif
return 0;
}
|
