aboutsummaryrefslogtreecommitdiffstats
path: root/test/ccured_olden/power/compute.c
diff options
context:
space:
mode:
Diffstat (limited to 'test/ccured_olden/power/compute.c')
-rw-r--r--test/ccured_olden/power/compute.c379
1 files changed, 0 insertions, 379 deletions
diff --git a/test/ccured_olden/power/compute.c b/test/ccured_olden/power/compute.c
deleted file mode 100644
index 83d47b59..00000000
--- a/test/ccured_olden/power/compute.c
+++ /dev/null
@@ -1,379 +0,0 @@
-/* For copyright information, see olden_v1.0/COPYRIGHT */
-
-/* compute.c
- *
- * By: Martin C. Carlisle
- * 6/15/94
- *
- * Implements computation phase of the Power Pricing problem
- * based on code by: Steve Lumetta, Sherry Li, and Ismail Khalil
- * University of California at Berkeley
- *
- */
-
-#include "power.h"
-#include <math.h>
-
-/*----------------------------------------------------------------------*/
-/* Leaf optimization 'global' variables */
-
- static double P=1.0;
- static double Q=1.0;
-
-
-/*----------------------------------------------------------------------*/
-/* Leaf optimization procedures */
-
-void optimize_node (double pi_R, double pi_I);
-double find_g ();
-double find_h ();
-double find_gradient_f (double pi_R, double pi_I, local double* gradient);
-double find_gradient_g (local double* gradient);
-double find_gradient_h (local double* gradient);
-void find_dd_grad_f (double pi_R, double pi_I, local double* dd_grad);
-double make_orthogonal (local double* v_mod, local double* v_static);
-
-
-void Compute_Tree(Root r) {
- register int i;
- Lateral l;
-#ifndef FUTURES
- Demand a;
-#else
- future_cell_demand fc[NUM_FEEDERS];
-#endif
- Demand tmp;
- double theta_R,theta_I;
-
- tmp.P = 0.0;
- tmp.Q = 0.0;
-#ifndef FUTURES
- for (i=0; i<NUM_FEEDERS; i++) {
- l = r->feeders[i];
- theta_R = r->theta_R;
- theta_I = r->theta_I;
- a = Compute_Lateral(l,theta_R,theta_I,theta_R,theta_I);
- tmp.P += a.P;
- tmp.Q += a.Q;
-
- }
-#else
- for (i=0; i<NUM_FEEDERS; i++) {
- l = r->feeders[i];
- theta_R = r->theta_R;
- theta_I = r->theta_I;
- FUTURE(l,theta_R,theta_I,theta_R,theta_I,Compute_Lateral,&fc[i]);
- }
- for (i=NUM_FEEDERS-1; i>=0; i--) {
- TOUCH(&fc[i]);
- tmp.P += fc[i].value.P;
- tmp.Q += fc[i].value.Q;
- }
-#endif
- r->D.P = tmp.P;
- r->D.Q = tmp.Q;
-}
-
-Demand Compute_Lateral(Lateral l, double theta_R, double theta_I,
- double pi_R, double pi_I) {
-#ifndef FUTURES
- Demand a1;
-#else
- future_cell_demand fc;
-#endif
- Demand a2;
- double new_pi_R, new_pi_I;
- double a,b,c,root;
- Lateral next;
- Branch br;
-
- new_pi_R = pi_R + l->alpha*(theta_R+(theta_I*l->X)/l->R);
- new_pi_I = pi_I + l->beta*(theta_I+(theta_R*l->R)/l->X);
-
- next = l->next_lateral;
- if (next != NULL)
-#ifndef FUTURES
- a1 = Compute_Lateral(next,theta_R,theta_I,new_pi_R,new_pi_I);
-#else
- FUTURE(next,theta_R,theta_I,new_pi_R,new_pi_I,Compute_Lateral,&fc);
-#endif
-
- br = l->branch;
- a2 = Compute_Branch(br,theta_R,theta_I,new_pi_R,new_pi_I);
-
- if (next != NULL) {
-#ifndef FUTURES
- l->D.P = a1.P + a2.P;
- l->D.Q = a1.Q + a2.Q;
-#else
- TOUCH(&fc);
- l->D.P = a2.P + fc.value.P;
- l->D.Q = a2.Q + fc.value.Q;
-#endif
- } else {
- l->D.P = a2.P;
- l->D.Q = a2.Q;
- }
-
- /* compute P,Q */
- a = l->R*l->R + l->X*l->X;
- b = 2*l->R*l->X*l->D.Q - 2*l->X*l->X*l->D.P - l->R;
- c = l->R*l->D.Q - l->X*l->D.P;
- c = c*c + l->R*l->D.P;
- root = (-b-sqrt(b*b-4*a*c))/(2*a);
- l->D.Q = l->D.Q + ((root-l->D.P)*l->X)/l->R;
- l->D.P = root;
-
- /* compute alpha, beta */
- a = 2*l->R*l->D.P;
- b = 2*l->X*l->D.Q;
- l->alpha = a/(1-a-b);
- l->beta = b/(1-a-b);
- return l->D;
-}
-
-Demand Compute_Branch(Branch br, double theta_R, double theta_I,
- double pi_R, double pi_I) {
- Demand a2,tmp;
- double new_pi_R, new_pi_I;
- double a,b,c,root;
- Leaf l;
- Branch next;
- int i;
-#ifdef FUTURES
- future_cell_demand fc;
-#else
- Demand a1;
-#endif
-
- new_pi_R = pi_R + br->alpha*(theta_R+(theta_I*br->X)/br->R);
- new_pi_I = pi_I + br->beta*(theta_I+(theta_R*br->R)/br->X);
-
- next = br->next_branch;
- if (next != NULL) {
-#ifndef FUTURES
- a1 = Compute_Branch(next,theta_R,theta_I,new_pi_R,new_pi_I);
-#else
- FUTURE(next,theta_R,theta_I,new_pi_R,new_pi_I,Compute_Branch,&fc);
-#endif
- }
-
- /* Initialize tmp */
- tmp.P = 0.0; tmp.Q = 0.0;
-
- for (i=0; i<LEAVES_PER_BRANCH; i++) {
- l = br->leaves[i];
- a2 = Compute_Leaf(l,new_pi_R,new_pi_I);
- tmp.P += a2.P;
- tmp.Q += a2.Q;
- }
- if (next != NULL) {
-#ifndef FUTURES
- br->D.P = a1.P + tmp.P;
- br->D.Q = a1.Q + tmp.Q;
-#else
- TOUCH(&fc);
- br->D.P = fc.value.P + tmp.P;
- br->D.Q = fc.value.Q + tmp.Q;
-#endif
- } else {
- br->D.P = tmp.P;
- br->D.Q = tmp.Q;
- }
-
- /* compute P,Q */
- a = br->R*br->R + br->X*br->X;
- b = 2*br->R*br->X*br->D.Q - 2*br->X*br->X*br->D.P - br->R;
- c = br->R*br->D.Q - br->X*br->D.P;
- c = c*c + br->R*br->D.P;
- root = (-b-sqrt(b*b-4*a*c))/(2*a);
- br->D.Q = br->D.Q + ((root-br->D.P)*br->X)/br->R;
- br->D.P = root;
- /* compute alpha, beta */
- a = 2*br->R*br->D.P;
- b = 2*br->X*br->D.Q;
- br->alpha = a/(1-a-b);
- br->beta = b/(1-a-b);
-
- return br->D;
-}
-
-Demand Compute_Leaf(Leaf l, double pi_R, double pi_I) {
- P = l->D.P;
- Q = l->D.Q;
-
- optimize_node(pi_R,pi_I);
-
- if (P<0.0) {
- P = 0.0;
- Q = 0.0;
- }
- l->D.P = P;
- l->D.Q = Q;
- return l->D;
-}
-
-/*----------------------------------------------------------------------*/
-
-void optimize_node (double pi_R, double pi_I)
-{
- double g;
- double h;
-
- double grad_f[2];
- double grad_g[2];
- double grad_h[2];
- double dd_grad_f[2];
- double magnitude;
-
- int i;
- double total;
- double max_dist;
-
- do {
- /* Move onto h=0 line */
- h=find_h ();
- if (fabs (h)>H_EPSILON) {
- magnitude=find_gradient_h (grad_h);
- total=h/magnitude;
- P-=total*grad_h[0];
- Q-=total*grad_h[1];
- }
-
- /* Check that g is still valid */
- g=find_g ();
- if (g>G_EPSILON) {
- magnitude=find_gradient_g (grad_g);
- find_gradient_h (grad_h);
- magnitude*=make_orthogonal (grad_g,grad_h);
- total=g/magnitude;
- P-=total*grad_g[0];
- Q-=total*grad_g[1];
- }
-
- /* Maximize benefit */
- magnitude=find_gradient_f (pi_R,pi_I,grad_f);
- find_dd_grad_f (pi_R,pi_I,dd_grad_f);
- total=0.0;
- for (i=0; i<2; i++)
- total+=grad_f[i]*dd_grad_f[i];
- magnitude/=fabs (total);
- find_gradient_h (grad_h);
- magnitude*=make_orthogonal (grad_f,grad_h);
- find_gradient_g (grad_g);
- total=0.0;
- for (i=0; i<2; i++)
- total+=grad_f[i]*grad_g[i];
- if (total>0) {
- max_dist=-find_g ()/total;
- if (magnitude>max_dist)
- magnitude=max_dist;
- }
- P+=magnitude*grad_f[0];
- Q+=magnitude*grad_f[1];
-
- h=find_h ();
- g=find_g ();
- find_gradient_f (pi_R,pi_I,grad_f);
- find_gradient_h (grad_h);
-
- } while (fabs (h)>H_EPSILON || g>G_EPSILON ||
- (fabs (g)>G_EPSILON &&
- fabs (grad_f[0]*grad_h[1]-grad_f[1]*grad_h[0])>F_EPSILON));
-}
-
-double find_g ()
-{
- return (P*P+Q*Q-0.8);
-}
-
-double find_h ()
-{
- return (P-5*Q);
-}
-
-double find_gradient_f (double pi_R, double pi_I, local double* gradient)
-{
- int i;
- double magnitude=0.0;
-
- gradient[0]=1/(1+P)-pi_R;
- gradient[1]=1/(1+Q)-pi_I;
- for (i=0; i<2; i++)
- magnitude+=gradient[i]*gradient[i];
- magnitude=sqrt (magnitude);
- for (i=0; i<2; i++)
- gradient[i]/=magnitude;
-
- return magnitude;
-}
-
-double find_gradient_g (local double* gradient)
-{
- int i;
- double magnitude=0.0;
-
- gradient[0]=2*P;
- gradient[1]=2*Q;
- for (i=0; i<2; i++)
- magnitude+=gradient[i]*gradient[i];
- magnitude=sqrt (magnitude);
- for (i=0; i<2; i++)
- gradient[i]/=magnitude;
-
- return magnitude;
-}
-
-double find_gradient_h (local double* gradient)
-{
- int i;
- double magnitude=0.0;
-
- gradient[0]=1.0;
- gradient[1]=-5.0;
- for (i=0; i<2; i++)
- magnitude+=gradient[i]*gradient[i];
- magnitude=sqrt (magnitude);
- for (i=0; i<2; i++)
- gradient[i]/=magnitude;
-
- return magnitude;
-}
-
-void find_dd_grad_f (double pi_R, double pi_I, local double* dd_grad)
-{
- double P_plus_1_inv=1/(P+1);
- double Q_plus_1_inv=1/(Q+1);
- double P_grad_term=P_plus_1_inv-pi_R;
- double Q_grad_term=Q_plus_1_inv-pi_I;
- double grad_mag;
-
- grad_mag=sqrt (P_grad_term*P_grad_term+Q_grad_term*Q_grad_term);
-
- dd_grad[0]=-P_plus_1_inv*P_plus_1_inv*P_grad_term/grad_mag;
- dd_grad[1]=-Q_plus_1_inv*Q_plus_1_inv*Q_grad_term/grad_mag;
-}
-
-double make_orthogonal (local double* v_mod, local double* v_static)
-{
- int i;
- double total=0.0;
- double length=0.0;
-
- for (i=0; i<2; i++)
- total+=v_mod[i]*v_static[i];
- for (i=0; i<2; i++) {
- v_mod[i]-=total*v_static[i];
- length+=v_mod[i]*v_mod[i];
- }
- length=sqrt (length);
- for (i=0; i<2; i++)
- v_mod[i]/=length;
-
- if (1-total*total<0) /* Roundoff error--vectors are parallel */
- return 0;
-
- return sqrt (1-total*total);
-}
-
-/*----------------------------------------------------------------------*/
generated by cgit (git 2.34.1) at 2024-06-29 04:36:05 +0000