Skip to content

Commit

Permalink
v479: better handling of V-S curves in lake-type reservoirs
Browse files Browse the repository at this point in the history 
-fixed how Q(h) and A(h) are recalculated for lake-type reservoirs with provided V(h) curves (ParseHRUFile.cpp/Reservoir.cpp)
-added CReservoir::GetReservoirName() (Reservoir.cpp/.h)
  • Loading branch information
James Craig authored and James Craig committed Sep 22, 2024
1 parent 1ce7fdd commit a6e2c7f
Show file tree
Hide file tree
Showing 4 changed files with 90 additions and 14 deletions.
45 changes: 38 additions & 7 deletions src/DemandOptimization.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -634,6 +634,15 @@ void CDemandOptimizer::InitializeDemands(CModel* pModel, const optStruct& Option

_demands_initialized=true;
}
string ComparisonToString(comparison C)
{
if (C==COMPARE_IS_EQUAL){return "="; }
if (C==COMPARE_GREATERTHAN){return ">"; }
if (C==COMPARE_LESSTHAN){return "<"; }
if (C==COMPARE_BETWEEN){return "in_range"; }
if (C==COMPARE_NOT_EQUAL){return "!="; }
return "?";
}
//////////////////////////////////////////////////////////////////
/// \brief Initializes Demand optimization instance
/// \notes to be called after .rvm file read
Expand Down Expand Up @@ -826,14 +835,14 @@ void CDemandOptimizer::InitializePostRVMRead(CModel* pModel, const optStruct& Op
cout<<" # Constraints/Goals: " <<_nGoals<<endl;
for (int i = 0; i < _nGoals; i++) {
if (_pGoals[i]->is_goal){tmpstr="[GOAL] "; }
else {tmpstr="[CONSTRAINT]"; }
else {tmpstr="[CONSTRAINT]"; }
cout<<" "<<i<<" "<<tmpstr<<": "<<_pGoals[i]->name<<endl;
for (int k=0; k<_pGoals[i]->nOperRegimes; k++)
{
comparison ctype=_pGoals[i]->pOperRegimes[k]->pExpression->compare;
cout<<" +oper regime: "<<_pGoals[i]->pOperRegimes[k]->reg_name<<endl;
cout<<" +expression: "<<_pGoals[i]->pOperRegimes[k]->pExpression->origexp<<endl;
cout<<" +comp. type: "<<ctype<<endl;
cout<<" +comp. type: "<<ComparisonToString(ctype)<<endl;
if (ctype==COMPARE_IS_EQUAL){
cout<<" +penalty: ("<<_pGoals[i]->pOperRegimes[k]->penalty_under<<","<<_pGoals[i]->pOperRegimes[k]->penalty_over<<")"<<endl;
}
Expand Down Expand Up @@ -1619,7 +1628,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio

Qin_last =pSB->GetReservoirInflow();
Qout_last=pRes->GetOutflowRate();
h_old =pRes->GetResStage();
h_old =pRes->GetResStage();//-h_ref
Adt =pRes->GetSurfaceArea()/Options.timestep/SEC_PER_DAY;
k =pRes->GetHRUIndex();
ET=0.0;
Expand Down Expand Up @@ -1669,9 +1678,9 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio
double Q_guess;
double h_sill;

Q_guess=pRes->GetDischargeFromStage(h_guess,nn);
dQdh =pRes->GetStageDischargeDerivative(h_guess,nn);
h_sill =pRes->GetSillElevation(nn);
Q_guess=pRes->GetDischargeFromStage(h_guess,nn);//+h_ref
dQdh =pRes->GetStageDischargeDerivative(h_guess,nn);//+h_ref
h_sill =pRes->GetSillElevation(nn);//-h_ref
if (_stage_discharge_as_goal)
{
// -----------------------------------------------------------------------
Expand Down Expand Up @@ -1718,7 +1727,7 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio
// Qout = 0 if h<h_sill
// handled via 6 constraint equations

const double LARGE_NUMBER=1e7;
const double LARGE_NUMBER=5000;

//------------------------------------------------------------------------
//Eqns 1 & 2 : h^n+1 - M*(1-b) <= h_sill
Expand Down Expand Up @@ -1938,14 +1947,34 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio
// Solve the LP problem!
// ----------------------------------------------------------------
//lp_lib::set_scaling(pLinProg, SCALE_CURTISREID+SCALE_EQUILIBRATE+SCALE_INTEGERS);
//lp_lib::set_scaling(pLinProg, SCALE_GEOMETRIC + SCALE_EQUILIBRATE + SCALE_INTEGERS +SCALE_DYNUPDATE);

//lp_lib::set_break_numeric_accuracy(pLinProg, 1e-6);
lp_lib::set_basiscrash(pLinProg, CRASH_MOSTFEASIBLE);

retval = lp_lib::solve(pLinProg);

// handle solve error
// ----------------------------------------------------------------
if (retval!=OPTIMAL)
{
cout<<"LP SOLVE CANNOT SOLVE OPTIMZIATION PROBLEM"<<endl;
cout<<"lp_lib::solve error code: "<<retval<<endl;
cout<<"date: "<<tt.date_string << endl;
cout<<"loop iteration: "<<iter<<" of "<<NUM_ITERATIONS<<endl;

/*lp_lib::get_variables (pLinProg,soln);
lp_lib::get_constraints(pLinProg,constr); //constraint matrix * soln
for (int j=0;j<nrows;j++)
{
_aSolverRowNames [j]=to_string(lp_lib::get_row_name(pLinProg,j+1));
_aSolverResiduals[j]=constr[j]-lp_lib::get_rh(pLinProg,j+1); // LHS-RHS for each goal/constraint
ctyp=lp_lib::get_constr_type(pLinProg,j+1);
if (ctyp==LE){upperswap(_aSolverResiduals[j],0.0); }
if (ctyp==GE){lowerswap(_aSolverResiduals[j],0.0); _aSolverResiduals[j]*=-1; }
cout<<j<<" "<<_aSolverRowNames [j]<<" : "<<_aSolverResiduals[j] << endl;
}*/
WriteLPSubMatrix(pLinProg,"overconstrained_lp_matrix.csv",Options);
ExitGracefully("SolveDemandProblem: non-optimal solution found. Problem is over-constrained. Remove or adjust management constraints.",RUNTIME_ERR);
}
Expand Down Expand Up @@ -1986,6 +2015,8 @@ void CDemandOptimizer::SolveDemandProblem(CModel *pModel, const optStruct &Optio
else if (typ == DV_QOUT ) {Q_iter[p]=value;}
}

//lp_lib::unscale(pLinProg);

// update lp matrix for iterative solution of stage-discharge curve and flow diversions
// ----------------------------------------------------------------
int hcol;
Expand Down
4 changes: 2 additions & 2 deletions src/ParseHRUFile.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1438,7 +1438,7 @@ CReservoir *ReservoirParse(CParser *p,string name,const CModel *pModel,long long
if(IsComment(s[0],Len)) { i--; }
else {
aV_ht[i] = s_to_d(s[0]);
aV[i] = s_to_d(s[1]);
aV [i] = s_to_d(s[1]);
}
}
p->Tokenize(s,Len); //:EndVolumeStageRelation
Expand Down Expand Up @@ -1924,7 +1924,7 @@ CReservoir *ReservoirParse(CParser *p,string name,const CModel *pModel,long long
}
if((type==CURVE_LAKE) && (aV!=NULL) && (aV_ht!=NULL))
{
pRes->SetVolumeStageCurve(aV_ht,aV,NV);//allows user to override prismatic lake assumption
pRes->SetVolumeStageCurve(aV_ht,aV,NV, weircoeff, cwidth);//allows user to override prismatic lake assumption
}
if((type==CURVE_LAKE) && (aA!=NULL) && (aA_ht!=NULL))
{
Expand Down
52 changes: 48 additions & 4 deletions src/Reservoir.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -398,6 +398,12 @@ CReservoir::~CReservoir()
//
long CReservoir::GetSubbasinID () const { return _SBID; }

//////////////////////////////////////////////////////////////////
/// \returns Reservoir name
//
string CReservoir::GetReservoirName() const { return _name; }


//////////////////////////////////////////////////////////////////
/// \returns reservoir storage [m3]
//
Expand Down Expand Up @@ -1003,16 +1009,54 @@ string CReservoir::GetCurrentConstraint() const {
/// \param nPoints - number of points in array
//

void CReservoir::SetVolumeStageCurve(const double *a_ht,const double *a_V,const int nPoints)
void CReservoir::SetVolumeStageCurve(const double *a_ht,const double *a_V,const int nPoints, double weircoeff, double crestwidth)
{
for(int i=0;i<_Np;i++)
for(int i=1;i<nPoints;i++)
{
_aVolume[i]=InterpolateCurve(_aStage[i],a_ht,a_V,nPoints,false);
if((i > 0) && ((_aVolume[i] - _aVolume[i-1]) <= -REAL_SMALL) && (_aVolume[i]!=0.0)) {
if (((a_ht[i]-a_ht[i-1])<=REAL_SMALL) || ((a_V[i]-a_V[i-1])<=REAL_SMALL)){
string warn = "CReservoir::SetVolumeStageCurve: volume-stage relationships must be monotonically increasing for all stages. [bad reservoir: " + _name + " "+to_string(_SBID)+"]";
ExitGracefully(warn.c_str(),BAD_DATA_WARN);
}
}
if (a_V[0]!=0.0){
string warn = "CReservoir::SetVolumeStageCurve: volume-stage relationships must start with volume of zero. [bad reservoir: " + _name + " "+to_string(_SBID)+"]";
ExitGracefully(warn.c_str(),BAD_DATA_WARN);
}

_min_stage =a_ht[0];
_max_stage =a_ht[nPoints-1];

//_Np=102; //SAME AS LAKE TYPE CONSTRUCTOR, BY NECESSITY!

//cout<<_name<<" CREST HEIGHT : "<<_crest_ht<<endl;
double dh;
string warn;
dh=(_max_stage-_min_stage)/(_Np-1);
_aStage [0]=_min_stage; //Memory
_aQ [0]=0.0;
_aQunder[0]=0.0;
_aArea [0]=0.0;
_aArea [0]=(InterpolateCurve(_aStage[0]+dh,a_ht,a_V,nPoints,false))/dh; //Area CANNOT go to zero
_aVolume[0]=0.0;
for (int i=1;i<_Np;i++)
{
_aStage [i]=_min_stage+i*dh;
_aQ [i]=0.0;
if ((_aStage[i]-_crest_ht)>0.0){
_aQ [i]=2.0/3.0*weircoeff*sqrt(2*GRAVITY)*crestwidth*pow((_aStage[i]-_crest_ht),1.5); //Overflow weir equation
}
if (((_aStage[i]-_crest_ht)<dh) && ((_aStage[i]-_crest_ht)>=0.0)){
_aStage[i]=_crest_ht; //ensures there is a point that perfectly coincides with crest
_aQ [i]=0.0;
}
_aQunder[i]=0.0;
_aVolume[i]=InterpolateCurve (_aStage[i] ,a_ht,a_V,nPoints,false);
_aArea [i]=(InterpolateCurve(_aStage[i]+dh,a_ht,a_V,nPoints,false)-_aVolume[i])/dh;
//cout<<" "<<_aStage[i]<<" "<<_aVolume[i]<<" "<<_aArea[i]<<" "<<_aQ [i]<<endl;

}
_max_capacity=_aVolume[_Np-1];

}
//////////////////////////////////////////////////////////////////
/// \brief overrides area stage curve for Lake-type reservoirs (if known)
Expand Down
3 changes: 2 additions & 1 deletion src/Reservoir.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -175,6 +175,7 @@ class CReservoir

//Accessors
long GetSubbasinID () const;
string GetReservoirName () const;

double GetStorage () const; //[m3]
double GetOldStorage () const; //[m3]
Expand Down Expand Up @@ -227,7 +228,7 @@ class CReservoir
void SetReservoirStage (const double &ht, const double &ht_last);
void SetControlFlow (const int i, const double &Q, const double &Qlast);
void SetOptimizedOutflow (const double &Qout);
void SetVolumeStageCurve (const double *a_ht,const double *a_V,const int nPoints);
void SetVolumeStageCurve (const double *a_ht,const double *a_V,const int nPoints, double weircoeff, double crestwidth);
void SetAreaStageCurve (const double *a_ht,const double *a_A,const int nPoints);
void SetGWParameters (const double &coeff, const double &h_ref);
void SetCrestWidth (const double &width);
Expand Down

0 comments on commit a6e2c7f

Please sign in to comment.