v512 - fixed demand optimization issues with lake at sill

Micromanaging managment support for at-sill lakes - A DEVIOUS BUG!
-Add CReservoir::GetSurfaceAreaFromStage() (Reservoir.h/cpp)
-Properly handle stage-discharge derivatives for h-h_sill<0.001 (Reservoir.cpp)
-(temporarily?) introduce DV_RSLACK for debugging (DemandOptimization.h; DemandExpressionHandling.cpp;DemandOptimization.cpp)
-track history of binary values through iterations (DemandOptimization.cpp)
-enable relaxation after first infeasible solution (DemandOptimization.cpp)
-fix bug in handling DSTAGE2 in iteration (DemandOptimization.cpp) [THIS DID THE TRICK]
-properly support changing lake area during iteration (DemandOptimization.cpp)
-don't update reservoir outflow unless managed at close of LP solve (DemandOptimization.cpp) [THE OTHER KEY TRICK]

minor:
Handle disabled SBs in CSubBasin::InitializePostRVM() (SubBasin.cpp)

Author: analytophile

src/CommonFunctions.cpp

@@ -1949,7 +1949,6 @@ double InterpolateCurve(const double x,const double *xx,const double *y,int N,bo
   {
     //int i=0; while ((x>xx[i+1]) && (i<(N-2))){i++;}//Dumb Search
     int i=SmartIntervalSearch(x,xx,N,ilast);
-    if(i==DOESNT_EXIST) { return 0.0; }
     ExitGracefullyIf(i==DOESNT_EXIST,"InterpolateCurve::mis-ordered list or infinite x",RUNTIME_ERR);
     ilast=i;
     if (fabs(xx[i+1]-xx[i]) < REAL_SMALL) { return (y[i]+y[i+1])/2; }  // x locations too close to each other

src/DemandExpressionHandling.cpp

@@ -288,7 +288,8 @@ string DVTypeToString(dv_type t)
   if (t==DV_QOUTRES ){return "DV_QOUTRES";  }
   if (t==DV_STAGE   ){return "DV_STAGE";    }
   if (t==DV_DSTAGE  ){return "DV_DSTAGE";   }
-  if (t==DV_DSTAGE2 ){return "DV_DSTAGE2";   }
+  if (t==DV_DSTAGE2 ){return "DV_DSTAGE2";  }
+  if (t==DV_RSLACK  ){return "DV_RSLACK";   }
   if (t==DV_BINRES  ){return "DV_BINRES";   }
   if (t==DV_DELIVERY){return "DV_DELIVERY"; }
   if (t==DV_RETURN  ){return "DV_RETURN";   }

src/DemandOptimization.cpp

@@ -44,6 +44,7 @@ enum dv_type
   DV_STAGE,   //< reservoir stage
   DV_DSTAGE,  //< positive change in reservoir stage over time step
   DV_DSTAGE2, //< negative change in reservoir stage over time step
+  DV_RSLACK,  //< slack variable for res outflow adjustment
   DV_BINRES,  //< binary integer value for above/beneath stage-discharge sill
   DV_DELIVERY,//< delivery of water demand
   DV_RETURN,  //< return flows to reach

src/DemandOptimization.h

@@ -1782,7 +1782,7 @@ bool ParseMainInputFile (CModel     *&pModel,
       break;
     }
     case(94):  //--------------------------------------------
-    {/*:AssimilatReservoirStage*/
+    {/*:AssimilateReservoirStage*/
       if(Options.noisy) { cout << "Assimilate lake stage on" << endl; }
       Options.assimilate_stage=true;
       break;

src/ParseInput.cpp

@@ -483,6 +483,11 @@ bool    CReservoir::UseInStageAssimilation() const
 //
 double  CReservoir::GetSurfaceArea       () const { return GetArea(_stage); }
 
+//////////////////////////////////////////////////////////////////
+/// \returns surface area given stage [m2]
+//
+double   CReservoir::GetSurfaceAreaFromStage  (const double h) const {return GetArea(h);}
+
 //////////////////////////////////////////////////////////////////
 /// \returns start-of-timestep surface area [m2]
 //
@@ -648,6 +653,8 @@ double CReservoir::GetStageDischargeDerivative(const double &stage, const int nn
   if(_pWeirHeightTS!=NULL) {
     weir_adj=_pWeirHeightTS->GetSampledValue(nn);
   }
+  double h_sill=GetSillElevation(nn);
+  if (stage <=h_sill){return 0.0;}
 
   return (GetWeirOutflow(stage+dh,weir_adj)-GetWeirOutflow(stage,weir_adj))/dh;
 }

src/Reservoir.cpp

@@ -193,6 +193,7 @@ class CReservoir
   double            GetOldStage              () const; //[m]
   double            GetSurfaceArea           () const; //[m2]
   double            GetOldSurfaceArea        () const; //[m2]
+  double            GetSurfaceAreaFromStage  (const double h) const;//[m2]
   double            GetLakebedThickness      () const; //[m]
   double            GetLakebedConductivity   () const; //[MJ/m/K/d]
   double            GetLakeConvectionCoeff   () const; //[MJ/m2/d/K]

src/Reservoir.h

@@ -1757,16 +1757,18 @@ void CSubBasin::Initialize(const double    &Qin_avg,          //[m3/s] from upst
 //
 void CSubBasin::InitializePostRVM(const optStruct &Options)
 {
-  if (_nWaterDemands>0){
-    _aQdelivered=new double [_nWaterDemands];
-    _aQreturned =new double [_nWaterDemands];
-    for (int ii=0;ii<_nWaterDemands;ii++){
-      _pWaterDemands[ii]->Initialize(Options);
-      _aQdelivered  [ii]=0.0;
-      _aQreturned   [ii]=0.0;
+  if (!_disabled) {
+    if (_nWaterDemands>0){
+      _aQdelivered=new double [_nWaterDemands];
+      _aQreturned =new double [_nWaterDemands];
+      for (int ii=0;ii<_nWaterDemands;ii++){
+        _pWaterDemands[ii]->Initialize(Options);
+        _aQdelivered  [ii]=0.0;
+        _aQreturned   [ii]=0.0;
+      }
     }
+    if (_pReservoir!=NULL){_pReservoir->InitializePostRVM(Options); }
   }
-  if (_pReservoir!=NULL){_pReservoir->InitializePostRVM(Options); }
 }
 //////////////////////////////////////////////////////////////////
 /// \brief Initializes channel flows, channel storage, rivulet storage from estimated mean inflows