/***************************************************************************
Synfuels Case Study
Formulation:   Continuous-Variable Approach
Max Function:  Grid Search or Golden-Section Search
February 7, 1997 by Jeff Stonebraker
filename:  SYNFUEL.C (compiled & run on MS QuickC Version 2.5)
***************************************************************************/

#include <stdio.h>
#include <math.h>
#include <float.h>

#define MAXBRANCH         6
#define MAXNODE          15
#define NUMSVARIABLES    14
#define PROGRAM           1
#define COST_1985         2
#define CARTEL_1985       3
#define FOREIGN_1985      4
#define EQUILIBRIUM_1985  5
#define EMBARGO_1985      6
#define EXPAND            7
#define COST_1995         8
#define CARTEL_1995       9
#define FOREIGN_1995     10
#define DEMAND_1995      11
#define EQUILIBRIUM_1995 12
#define EMBARGO_1995     13
#define COST_FACTOR      14
#define CHANCE            1
#define DECISION          2
#define ENDPOINT          3
#define AUXILIARY         4
#define INFINITY          1.0E30
#define r                 (sqrt(5.0)-1.0)/2.0
#define MAXGRID           0
#define MAXGOLDEN         1

double EU(int i, double s[NUMSVARIABLES+1]);
int n(int i, int j, double s[NUMSVARIABLES+1]);
double p(int i, int j, double s[NUMSVARIABLES+1]);
double f(int i, int j, double s[NUMSVARIABLES+1]);
double u(double x);
double CF(double x, double y);
double maxgrid(int i, double s[NUMSVARIABLES+1]);
double maxgolden(int i, double s[NUMSVARIABLES+1]);
double L(int i, double s[NUMSVARIABLES+1]);
double U(int i, double s[NUMSVARIABLES+1]);
double Xmax[NUMSVARIABLES+1];
double s[NUMSVARIABLES+1];
double lo[NUMSVARIABLES+1];
double hi[NUMSVARIABLES+1];
double accuracy[NUMSVARIABLES+1];
double N(double mu, double sigma, int j);
int b[MAXNODE+1];
int t[MAXNODE+1];
int v[MAXNODE+1];
FILE *fp_out;

main()
{
   double EUmax;
   b[2] =  3;
   b[3] =  2;
   b[4] =  3; 
   b[6] =  2; 
   b[8] =  3;
   b[9] =  2; 
   b[10] = 3; 
   b[11] = 3; 
   b[13] = 2;
   t[1] =  DECISION; 
   t[2] =  CHANCE; 
   t[3] =  CHANCE; 
   t[4] =  CHANCE;
   t[5] =  AUXILIARY; 
   t[6] =  CHANCE; 
   t[7] =  DECISION;
   t[8] =  CHANCE; 
   t[9] =  CHANCE; 
   t[10] = CHANCE; 
   t[11] = CHANCE;
   t[12] = AUXILIARY; 
   t[13] = CHANCE; 
   t[14] = AUXILIARY; 
   t[15] = ENDPOINT;
   v[1] =  PROGRAM; 
   v[2] =  COST_1985; 
   v[3] =  CARTEL_1985; 
   v[4] =  FOREIGN_1985;
   v[5] =  EQUILIBRIUM_1985; 
   v[6] =  EMBARGO_1985; 
   v[7] =  EXPAND; 
   v[8] =  COST_1995;
   v[9] =  CARTEL_1995; 
   v[10] = FOREIGN_1995; 
   v[11] = DEMAND_1995;
   v[12] = EQUILIBRIUM_1995; 
   v[13] = EMBARGO_1995; 
   v[14] = COST_FACTOR;
   accuracy[v[1]] = 0.586/4.0;
   accuracy[v[7]] = 1.825/4.0;
   fp_out = fopen("synfuel.out", "a");
   #if MAXGRID == 1
      EUmax = maxgrid(1, s);
      fprintf(fp_out, "\n\nPROGRAM = %.5f and EV = %.5f", Xmax[v[1]], EUmax);
      fprintf(fp_out, "\n\tGrid Search");
      fprintf(fp_out, "\n\tPROGRAM accuracy = %.5f & EXPAND accuracy = %.5f",
		accuracy[v[1]], accuracy[v[7]]);
   #endif
   #if MAXGOLDEN == 1
      EUmax = maxgolden(1, s);
      fprintf(fp_out, "\n\nPROGRAM = %.5f and EV = %.5f", Xmax[v[1]], EUmax);
      fprintf(fp_out, "\n\tGolden-Section Search");
      fprintf(fp_out, "\n\tPROGRAM accuracy = %.5f & EXPAND accuracy = %.5f",
		accuracy[v[1]], accuracy[v[7]]);
   #endif
   fclose(fp_out);
}

/*************************************************************************** 
Expected Utility (EU) Function for chance node (expectation operator), 
decision node (maximization operator), endpoint node (utility operator), and
auxiliary node (deterministic function). 
***************************************************************************/

double EU(int i, double s[NUMSVARIABLES+1])
{
   int j;
   double NewEU;
   switch(t[i])
   {
      case CHANCE:          
	 NewEU = 0.0;
	 for(j = 1; j <= b[i]; ++j)
	 {
	    s[v[i]] = f(i, j, s);
	    NewEU += p(i, j, s)*EU(n(i, j, s), s);
	 }
	 return(NewEU);
	 break;             
      case DECISION:
	 #if MAXGRID == 1
	    NewEU = maxgrid(i, s);
	 #endif
	 #if MAXGOLDEN == 1
	    NewEU = maxgolden(i, s);
	 #endif
	 return(NewEU);
	 break;
      case ENDPOINT:        
	 return(u(f(i, 1, s)));
	 break;
      case AUXILIARY:       
	 s[v[i]] = f(i, 1, s);
	 return(EU(n(i, 1, s), s));
	 break;
   }
}

/*************************************************************************** 
Maximization procedure for a decision node i using discrete-grid search
(maxgrid).  Initializes lower (lo[v[i]]) and upper (hi[v[i]]) bound values
for decision node i.  Determines the number of discrete increments for
decison node i (dummy_b) using the lower & upper bound values and the
desired accuracy level of decision node i (accuracy[v[i]]).  Generates 
discrete increments for decision node i (s[v[i]), beginning at its lower 
bound and continuing to its upper bound in increments of the desired 
accuracy level.  Computes the expected utility of each next node for each 
discrete increment of decision node i (testEU). Determines and records 
optimal solution (EUmax and Xmax[v[i]].)
***************************************************************************/

double maxgrid(int i, double s[NUMSVARIABLES+1])
{
   double testEU, EUmax, dummy_b; int j;
   EUmax = -INFINITY;
   lo[v[i]] = L(i, s);
   hi[v[i]] = U(i, s); 
   dummy_b = ceil((hi[v[i]]-lo[v[i]])/accuracy[v[i]])+1.0;
   b[i] = (int)(dummy_b);
   for(j = 1; j < b[i]; ++j)
   {
      s[v[i]] = (double)(j-1)*accuracy[v[i]]+lo[v[i]];
      testEU = EU(n(i, j, s), s);
      if(testEU > EUmax)
      {
	 EUmax = testEU;
	 Xmax[v[i]] =  s[v[i]];
      }
   }
   return(EUmax);
}

/***************************************************************************
Maximization procedure for a decision node i using the golden-section search
(maxgolden). Initializes lower (lo[v[i]]) and upper (hi[v[i]]) bound values 
for decision node i.  The initial if statement tests whether the 
golden-section search is done.  If it is done, then it updates the value 
(Xmax[v[i]]) and expected utility (EUmax) of the optimal solution else it
determines the total number of search intervals (iteration).  Determines 
trial solutions (x and y) using lower & upper bound values of decision node 
i, and the golden-section ratio (r).  Then the search interval narrows until 
within the desired accuracy level of decision node i (accuracy[v[i]]).  
Finally, determines and records the value and expected utility of the 
optimal solution.
***************************************************************************/

double maxgolden(int i, double s[NUMSVARIABLES+1])
{
   double x, y, testEU_lo, testEU_hi, EUmax;
   int j, k, iteration;
   lo[v[i]] = L(i, s);
   hi[v[i]] = U(i, s);
   if((hi[v[i]] - lo[v[i]]) == 0.0)
   {
      Xmax[v[i]] = (lo[v[i]]+hi[v[i]])/2.0;
      s[v[i]] = Xmax[v[i]];
      EUmax = EU(n(i, j, s), s);
      return(EUmax);
   }
   else
   {
      iteration = (int)(ceil(log(accuracy[v[i]]/(hi[v[i]]-lo[v[i]]))/log(r)));
      x = (hi[v[i]]-lo[v[i]])*r+lo[v[i]];
      y = (hi[v[i]]-lo[v[i]])*(1.0-r)+lo[v[i]];
      s[v[i]] = x;
      testEU_hi = EU(n(i, j, s), s);
      s[v[i]] = y;
      testEU_lo = EU(n(i, j, s), s);
      for(k = 1; k <= iteration - 1; ++k)
      {
	 if(testEU_hi >= testEU_lo)
	 {
	    lo[v[i]] = y;
	    y = x;
	    x = (hi[v[i]]-lo[v[i]])*r+lo[v[i]];
	    s[v[i]] = x;
	    testEU_lo = testEU_hi;
	    testEU_hi = EU(n(i, j, s), s);
	 }
	 else
	 {
	    hi[v[i]] = x;
	    x = y;
	    y = (hi[v[i]]-lo[v[i]])*(1.0-r)+lo[v[i]];
	    s[v[i]] = y;
	    testEU_hi = testEU_lo;
	    testEU_lo = EU(n(i, j, s), s);
	 }
      }
      Xmax[v[i]] = (lo[v[i]]+hi[v[i]])/2.0;
      s[v[i]] = Xmax[v[i]];
      EUmax = EU(n(i, j, s), s);
      return(EUmax);
   }
}

/***************************************************************************
Upper Bound Function (U) sets a decision node i to its upper bound value.
***************************************************************************/

double U(int i, double s[NUMSVARIABLES+1])
{
   switch(i)
   {
      case 1:
	 hi[v[1]] = 0.586;
	 break;
      case 7:
	 hi[v[7]] = 1.825;
	 break;
   }
   return(hi[v[i]]);
}

/***************************************************************************
Lower Bound Function (L) sets a decision node i to its lower bound value.
***************************************************************************/

double L(int i, double s[NUMSVARIABLES+1])
{
   lo[v[i]] = 0.0;
   return(lo[v[i]]);
}

/***************************************************************************
The Next Node Function (n) provides the algebraic structure of the
sequential decision model (decision tree).
***************************************************************************/

int n(int i, int j, double s[NUMSVARIABLES+1])
{
   return(i+1);
}

/***************************************************************************
The Utility Function (u) can be used to introduce risk attitude.  In this
example, the risk attitude is risk neutral or expected value.
***************************************************************************/

double u(double x)
{
   return(x);
}

/***************************************************************************
Branch Probability Functions (p) are used to represent the probability of a
chance node in a decision model.
***************************************************************************/

double p(int i, int j, double s[NUMSVARIABLES+1])
{
   switch(i)
   {
/***************************************************************************
Cases 2, 4, 8, 10, and 11 are Branch Probability Functions for Cost 1985,
Foreign 1985, Cost 1995, Foreign 1995, and Demand 1995, respectively.  They
are also assumed to be continuous random variables and normally distributed.  
These variables were integrated using the extended Pearson-Tukey (EP-T) 
approximation.  The probabilities for the 0.05, 0.50, and 0.95 fractile
values of the random variable, using the EP-T approximation, are 0.185, 
0.630, and 0.185, respectively.
***************************************************************************/
      case 2:     
      case 4:     
      case 8:     
      case 10:    
      case 11:     
	 switch(j)
	 {
	    case 1:
	       return(0.185);   
	       break;           
	    case 2:
	       return(0.630);  
	       break;          
	    case 3:
	       return(0.185);  
	       break;          
	 }
	 break;
/***************************************************************************
Case 3 is the Branch Probability Function for Cartel 1985.
***************************************************************************/
      case 3:      
	 return(0.5);
	 break;
/***************************************************************************
Case 6 is the Branch Probability Function for Embargo 1985.
***************************************************************************/
      case 6:    
	 switch(j)
	 {
	    case 1:
	       return(0.1);   
	       break;
	    case 2:
	       return(0.9);   
	       break;
	 }
	 break;
/***************************************************************************
Case 9 is the Branch Probability Function for Cartel 1995 given the outcome
of the cartel in 1985 (case 1 is strong and case 2 is weak).
***************************************************************************/
      case 9:     
	 switch((int)(s[CARTEL_1985]))
	 {
	    case 1:                    
	       switch(j)
	       {
		  case 1:              
		     return(0.8);
		     break;
		  case 2:              
		     return(0.2);
		     break;
	       }
	       break;
	    case 2:                     
	       switch(j)
	       {
		  case 1:              
		     return(0.2);
		     break;
		  case 2:              
		     return(0.8);
		     break;
	       }
	       break;
	 }
	 break;
/***************************************************************************
Case 13 is the Branch Probability Function for Embargo 1995.
***************************************************************************/
      case 13:     
	 switch(j)
	 {
	    case 1:             
	       return(0.05);
	       break;
	    case 2:
	       return(0.95);   
	       break;
	 }
	 break;
   }
}

/***************************************************************************
Branch Value Functions (f) are used to represent the value of a chance,
auxiliary, or endpoint node in a decision model.
***************************************************************************/

double f(int i, int j, double s[NUMSVARIABLES+1])
{
   double market_clearing_price, c, mu, sigma, fractile_value;
   switch(i)
   {
/***************************************************************************
Case 2 is the Branch Value Function for Cost 1985.
***************************************************************************/
      case 2:
	 mu = 17.7; sigma = 3.74;
	 fractile_value = N(mu, sigma, j);
	 return (fractile_value);
	 break;
/***************************************************************************
Cases 3 and 9 are Branch Value Functions for Cartel 1985 and Cartel 1995.
***************************************************************************/
      case 3: case 9:
	 return((double)(j));
	 break;
/***************************************************************************
Case 4 is the Branch Value Function for Foreign 1985.
***************************************************************************/
      case 4:           
	 switch((int)(s[CARTEL_1985]))
	 {
	    case 1:
	       mu = 15.0; sigma = 2.83;
	       fractile_value = N(mu, sigma, j);
	       return (fractile_value);
	       break;
	    case 2:
	       mu = 8.0; sigma = 1.41;
	       fractile_value = N(mu, sigma, j);
	       return (fractile_value);
	       break;
	 }
	 break;
/***************************************************************************
Case 5 is the Branch Value Function for Equilibrium 1985.
***************************************************************************/
      case 5:         
	 market_clearing_price = (1888.875/(19.893+s[PROGRAM]))-69.0;
	 if(market_clearing_price <= s[FOREIGN_1985])
	 {
	    return(market_clearing_price);
	 }
	 else
	 {
	    return(s[FOREIGN_1985]);
	 }
	 break;
/***************************************************************************
Case 6 is the Branch Value Function for Embargo 1985.
***************************************************************************/
      case 6:          
	 switch(j)
	 {
	    case 1:
	       c = (1888.875/(s[EQUILIBRIUM_1985]+69.0))-19.893;
	       return(-1967.58*(c*c-(0.5*(c+s[PROGRAM]))*(0.5*(c+s[PROGRAM])))/
		      ((19.893+c)*(19.893+c)));
	       break;
	    case 2:
	       return(0.0);
	       break;
	 }
	 break;
/***************************************************************************
Case 8 is the Branch Value Function for Cost 1995.
***************************************************************************/
      case 8:           
	 mu = 16.9 - 7.54*s[PROGRAM];
	 sigma = 3.56 - 2.14*s[PROGRAM];
	 fractile_value = N(mu, sigma, j);
	 return (fractile_value);
	 break;
/***************************************************************************
Case 10 is the Branch Value Function for Foreign 1995.
***************************************************************************/
      case 10:        
	 switch((int)(s[CARTEL_1995]))
	 {
	    case 1:
	       if((s[COST_1995] >= 5.0) && (s[COST_1995] <= 14.0))
	       {
		  mu = 9.39+0.542*s[COST_1995];
		  sigma = 2.36+0.135*s[COST_1995];
	       }
	       else
	       {
		  mu = 14.6+0.175*s[COST_1995];
		  sigma = 3.41+0.0590*s[COST_1995];
	       }
	       fractile_value = N(mu, sigma, j);
	       return (fractile_value);
	       break;
	    case 2:
	       mu = 11.3;
	       sigma = 2.49;
	       fractile_value = N(mu, sigma, j);
	       return (fractile_value);
	       break;
	 }
	 break;
/***************************************************************************
Case 11 is the Branch Value Function for Demand 1995.
***************************************************************************/
      case 11:  
	 mu = 15.4;
	 sigma = 2.07;
	 fractile_value = N(mu, sigma, j);
	 return (fractile_value);
	 break;
/***************************************************************************
Case 12 is the Branch Value Function for Equilibrium 1995.
***************************************************************************/
      case 12:      
	 market_clearing_price = (809.910/(s[DEMAND_1995]
	 +s[PROGRAM]+s[EXPAND]))-23.96;
	 if(market_clearing_price <= s[FOREIGN_1995])
	 {
	    return(market_clearing_price);
	 }
	 else
	 {
	    return(s[FOREIGN_1995]);
	 }
	 break;        
/***************************************************************************
Case 13 is the Branch Value Function for Embargo 1995.
***************************************************************************/
      case 13:          
	 switch(j)
	 {
	    case 1:
	       c = (809.910/(s[EQUILIBRIUM_1995]+23.946))-s[DEMAND_1995];
	       return(-843.66*(c*c-(0.5*(c+s[PROGRAM]+s[EXPAND]))*
		      (0.5*(c+s[PROGRAM]+s[EXPAND])))/
		      ((s[DEMAND_1995]+c)*(s[DEMAND_1995]+c)));
	       break;
	    case 2:
	       return(0.0);
	       break;
	 }
	 break;
/***************************************************************************
Case 14 is the Branch Value Function for Cost Factor.
***************************************************************************/
      case 14:          
	 if(s[EXPAND] >= 0.0 && s[EXPAND] <= 0.365)
	 {
	    return(CF(s[PROGRAM], 0.000)*(s[PROGRAM]+s[EXPAND]));
	 }
	 else if(s[EXPAND] >= 0.365 && s[EXPAND] <= 0.730)
	 {
	    return(0.1825*(CF(s[PROGRAM], 0.365)+CF(s[PROGRAM], 0.000))
	    +CF(s[PROGRAM], 0.365)*(s[PROGRAM]+s[EXPAND]-0.365));
	 }
	 else if(s[EXPAND] >= 0.730 && s[EXPAND] <= 1.095)
	 {
	    return(0.1825*(CF(s[PROGRAM], 0.000)+2.0*CF(s[PROGRAM], 0.365)
	    +CF(s[PROGRAM], 0.730))
	    +CF(s[PROGRAM], 0.730)*(s[PROGRAM]+s[EXPAND]-0.730));
	 }
	 else if(s[EXPAND] >= 1.095 && s[EXPAND] <= 1.460)
	 {
	    return(0.1825*(CF(s[PROGRAM], 0.000)+2.0*CF(s[PROGRAM], 0.365)
	    +2.0*CF(s[PROGRAM], 0.730)+CF(s[PROGRAM], 1.095))
	    +CF(s[PROGRAM], 1.095)*(s[PROGRAM]+s[EXPAND]-1.095));
	 }
	 else
	 {
	    return(0.1825*(CF(s[PROGRAM], 0.000)+2.0*CF(s[PROGRAM], 0.365)
	    +2.0*CF(s[PROGRAM], 0.730)+2.0*CF(s[PROGRAM], 1.095)
	    +CF(s[PROGRAM], 1.460))
	    +CF(s[PROGRAM], 1.460)*(s[PROGRAM]+s[EXPAND]-1.460));
	 }
	 break;
/***************************************************************************
Case 15 is the Branch Value (Endpoint) Function for Net Benefit.
***************************************************************************/
      case 15:          
	 return(4.2*(
		     8084.49-1888.875*log(s[EQUILIBRIUM_1985]+69.0)
		     +19.893*s[EQUILIBRIUM_1985]+s[EMBARGO_1985]
		     +(s[EQUILIBRIUM_1985]-s[COST_1985])*s[PROGRAM]
		     -0.4*s[PROGRAM]
		     )+
		1.62*(
		     809.91*(log(809.91/s[DEMAND_1995])-1.0)
		     -809.91*log(s[EQUILIBRIUM_1995]+23.946)
		     +s[DEMAND_1995]*(s[EQUILIBRIUM_1995]+23.946)
		     +s[EMBARGO_1995]+s[EQUILIBRIUM_1995]*(s[PROGRAM]+s[EXPAND])
		     -s[COST_1995]*s[COST_FACTOR]
		     +(s[EQUILIBRIUM_1995]-s[COST_1985])*s[PROGRAM]
		     -0.4*(s[PROGRAM]+s[EXPAND])
		     )
	       );
	 break;
   }
}

/***************************************************************************
Cost 1985, Cost 1995, Foreign 1985, Foreign 1995, and Demand 1995 are 
assumed to be continuous random variables and normally distributed (mu and
sigma).  This function (N) returns the branch values for the 0.05, 0.50, and 
0.95 fractiles of these variables using extended Pearson-Tukey (EP-T) 
approximation.
***************************************************************************/

double N(double mu, double sigma, int j)
{         
    switch(j)
	{
	    case 1:
	       return(mu - 1.645*sigma);  
	       break;                     
	    case 2:
	       return(mu);                
	       break;                     
	    case 3:
	       return(mu + 1.645*sigma);  
	       break;                     
	 }
}

/***************************************************************************
This function (CF) linear interpolates values between adjacent Program and 
Expand entries for the capacity expansion cost factors. 
***************************************************************************/

double CF(double x, double y)
{
   if(x >=  0.0 && x <= 0.115)
   {
      if(y == 0.000)
      {
	 return(0.94-0.348*x);
      }
      else if(y == 0.365)
      {
	 return(1.0-0.348*x);
      }
      else if(y == 0.730)
      {
	 return(1.13-0.696*x);
      }
      else if(y == 1.095)
      {
	 return(1.32-1.39*x);
      }
      else if(y == 1.460)
      {
	 return(3.0-14.7*x);
      }
      else
      {
	 return(10.0-73.0*x);
      }
   }
   else if(x >= 0.115 && x <= 0.339)
   {
      if(y == 0.000)
      {
	 return(0.905-0.0446*x);
      }
      else if(y == 0.365)
      {
	 return(0.97-0.0893*x);
      }
      else if(y == 0.730)
      {
	 return(1.07-0.179*x);
      }
      else if(y == 1.095)
      {
	 return(1.19-0.223*x);
      }
      else if(y == 1.460)
      {
	 return(1.34-0.313*x);
      }
      else
      {
	 return(1.69-0.67*x);
      }
   }
   else
   {
      if(y == 0.000)
      {
	 return(0.904-0.0405*x);
      }
      else if(y == 0.365)
      {
	 return(0.954-0.0405*x);
      }
      else if(y == 0.730)
      {
	 return(1.04-0.081*x);
      }
      else if(y == 1.095)
      {
	 return(1.18-0.202*x);
      }
      else if(y == 1.460)
      {
	 return(1.33-0.283*x);
      }
      else
      {
	 return(1.6-0.405*x);
      }
   }
}