@@ -357,7 +357,7 @@ void CDemandOptimizer::AddUserDecisionVar(const decision_var* pDV)
357357void CDemandOptimizer::OverrideSDCurve (const int p)
358358{
359359 _aDisableSDCurve[p]=true ;
360- }
360+ }
361361// ////////////////////////////////////////////////////////////////
362362// / \brief sets bounds for user specified decision variable
363363//
@@ -1375,6 +1375,22 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
13751375 int *lpsbrow=new int [_pModel->GetNumSubBasins ()]; // index of constraint equation for subbasin reaches
13761376 int *lpgoalrow=new int [_nGoals]; // index of goal equation for all user-specified goals
13771377
1378+ double **aQiterHist=NULL ,**ahiterHist=NULL ;
1379+
1380+ aQiterHist=new double *[_pModel->GetNumSubBasins ()];
1381+ ahiterHist=new double *[_pModel->GetNumSubBasins ()];
1382+ ExitGracefullyIf (ahiterHist==NULL ," CDemandOptimizer::InitializePostRVMRead" ,OUT_OF_MEMORY);
1383+ for (int p=0 ;p<_pModel->GetNumSubBasins ();p++){
1384+ ahiterHist[p]=NULL ;
1385+ aQiterHist[p] = new double [_maxIterations];
1386+ ahiterHist[p] = new double [_maxIterations];
1387+ ExitGracefullyIf (ahiterHist[p]==NULL ," CDemandOptimizer::InitializePostRVMRead (2)" ,OUT_OF_MEMORY);
1388+ for (int i = 0 ; i < _maxIterations; i++) {
1389+ aQiterHist[p][i]=0.0 ;
1390+ ahiterHist[p][i]=0.0 ;
1391+ }
1392+ }
1393+
13781394 // instantiate linear programming solver
13791395 // ----------------------------------------------------------------
13801396 lp_lib::lprec *pLinProg;
@@ -1807,7 +1823,6 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
18071823
18081824 IncrementAndSetRowName (pLinProg,rowcount," dh_def_" +to_string (pSB->GetID ()));
18091825
1810-
18111826 if (!_aDisableSDCurve[p])
18121827 {
18131828 // Second goal equation (relation between Qout and h):
@@ -1852,6 +1867,12 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
18521867 ExitGracefullyIf (retval==0 ," SolveManagementProblem::Error adding stage discharge constraint B" ,RUNTIME_ERR);
18531868 IncrementAndSetRowName (pLinProg,rowcount," reserv_Q_B" +to_string (pSB->GetID ()));
18541869
1870+ // ------------------------------------------------------------------------
1871+ // Eqn 3 : Q^n+1 + M*(b) <= M
1872+ // M must be >> Q^n+1
1873+ // enforces Q^n+1==0 when b==1
1874+ // ------------------------------------------------------------------------
1875+
18551876 col_ind[0 ]=GetDVColumnInd (DV_QOUTRES,_aResIndices[p]); row_val[0 ]=1.0 /sqrt (LARGE_NUMBER);
18561877 col_ind[1 ]=GetDVColumnInd (DV_BINRES ,_aResIndices[p]); row_val[1 ]=+sqrt (LARGE_NUMBER);
18571878 RHS =+sqrt (LARGE_NUMBER);
@@ -1861,10 +1882,10 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
18611882 IncrementAndSetRowName (pLinProg,rowcount," reserv_Q_D" +to_string (pSB->GetID ()));
18621883
18631884 // ------------------------------------------------------------------------
1864- // Eqns 5 & 6 : Q^n+1 >= Qout_guess + dQ/dh * (h^n+1-h_guess) - M(b)
1885+ // Eqns 4 & 5 : Q^n+1 >= Qout_guess + dQ/dh * (h^n+1-h_guess) - M(b)
18651886 // Q^n+1 <= Qout_guess + dQ/dh * (h^n+1-h_guess) + M(b)
1866- // Q^n+1 - dQdh * h^n+1 + Mb >= [Qout_guess - dQ/dh * h_guess]
1867- // Q^n+1 - dQdh * h^n+1 - Mb <= [Qout_guess - dQ/dh * h_guess]
1887+ // Q^n+1 - dQ/dh * h^n+1 + Mb >= [Qout_guess - dQ/dh * h_guess]
1888+ // Q^n+1 - dQ/dh * h^n+1 - Mb <= [Qout_guess - dQ/dh * h_guess]
18681889 // dQ << M
18691890 // ------------------------------------------------------------------------
18701891
@@ -1889,6 +1910,8 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
18891910 ExitGracefullyIf (retval==0 ," SolveManagementProblem::Error adding stage discharge constraint F" ,RUNTIME_ERR);
18901911 IncrementAndSetRowName (pLinProg,rowcount," reserv_Q_F" +to_string (pSB->GetID ()));
18911912
1913+ Q_iter[p]=Q_guess;
1914+ h_iter[p]=h_guess;
18921915 }
18931916 else /* if (aDisableSDCurve[p])*/ // keep same rows - make inert equations
18941917 {
@@ -1905,6 +1928,9 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
19051928 retval = lp_lib::add_constraintex (pLinProg,1 ,row_val,col_ind,ROWTYPE_LE,RHS);
19061929 IncrementAndSetRowName (pLinProg,rowcount," reserv_Q_F" +to_string (pSB->GetID ()));
19071930
1931+ Q_iter[p]=-1 ;
1932+ h_iter[p]=-1 ;
1933+
19081934 s+=2 ; // to ensure EnvMin counter is working
19091935 }
19101936 res_count++;
@@ -2024,11 +2050,6 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
20242050 double norm;
20252051 do
20262052 {
2027- if (_do_debug_level==2 )// EXTREME OUTPUT!!
2028- {
2029- WriteLPSubMatrix (pLinProg," lp_matrix.csv" ,Options);
2030- }
2031-
20322053 // ----------------------------------------------------------------
20332054 // Solve the LP problem!
20342055 // ----------------------------------------------------------------
@@ -2038,6 +2059,7 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
20382059 // lp_lib::set_break_numeric_accuracy(pLinProg, 1e-6);
20392060 // lp_lib::set_basiscrash(pLinProg, CRASH_MOSTFEASIBLE);
20402061 // lp_lib::set_verbose(pLinProg, DETAILED); //FOR DEBUGGING - BUT THIS IS WAY WAY TOO MUCH INFO
2062+ // lp_lib::set_epspivot(pLinProg, 1e-9);
20412063 retval = lp_lib::solve (pLinProg);
20422064 if (retval!=OPTIMAL)
20432065 {
@@ -2084,34 +2106,62 @@ void CDemandOptimizer::SolveManagementProblem(CModel *pModel, const optStruct &O
20842106 cout<<" lp_lib::solve error code: " <<retval<<" (" <<code<<" )" <<endl;
20852107 cout<<" lp accuracy: " <<lp_lib::get_accuracy (pLinProg)<<endl;
20862108 cout<<" date: " <<tt.date_string << endl;
2087- cout<<" loop iteration: " <<iter<<" of " <<_maxIterations<<endl;
2109+ cout<<" loop iteration: " <<iter+ 1 <<" of " <<_maxIterations<<endl;
20882110
20892111 cout<<" Problematic Constraints? (largest residuals): " <<endl;
20902112 for (int j=0 ;j<nrows;j++)
20912113 {
20922114 if (fabs (_aSolverResiduals[j])>REAL_SMALL)
2093- cout<<" - " <<_aSolverRowNames[j] << " " << _aSolverResiduals[j] << endl;
2115+ cout<<" > " <<_aSolverRowNames[j] << " " << _aSolverResiduals[j] << endl;
20942116 }
20952117
2096- // lp_lib::print_debugdump(pLinProg, dumpfile.c_str());
2118+ cout<<" Iteration Sequence: " <<endl;
2119+ for (int p = 0 ; p < pModel->GetNumSubBasins (); p++) {
2120+ if (pModel->GetSubBasin (p)->GetReservoir ()!=NULL ){
2121+ long long SBID=pModel->GetSubBasin (p)->GetID ();
2122+ ahiterHist[p][iter]=h_iter[p];
2123+ aQiterHist[p][iter]=Q_iter[p];
2124+ cout << " >Q " << SBID << " : " ; for (int i = 0 ; i <= iter; i++) { cout << aQiterHist[p][i] << " " ; }cout << endl;
2125+ cout << " >h " << SBID << " : " ; for (int i = 0 ; i <= iter; i++) { cout << ahiterHist[p][i] << " " ; }cout << endl;
2126+ }
2127+ }
2128+ char dump[256 ];
2129+ dump[0 ] = ' \0 ' ;
2130+ strcpy (dump, dumpfile.c_str ());
2131+ lp_lib::print_debugdump (pLinProg,dump);
20972132 WriteLPSubMatrix (pLinProg," overconstrained_lp_matrix.csv" ,Options);
2133+
2134+ lp_lib::delete_lp (pLinProg);
20982135 ExitGracefully (" SolveManagementProblem: non-optimal solution found. Problem is over-constrained. Remove or adjust management constraints." ,RUNTIME_ERR);
20992136 }
2137+ if (_do_debug_level==2 )// EXTREME OUTPUT!! - last valid lp
2138+ {
2139+ WriteLPSubMatrix (pLinProg," lp_matrix.csv" ,Options);
2140+ }
21002141
21012142 // grab and store improved estimate of reservoir stage + flow diversion
21022143 // ----------------------------------------------------------------
21032144 double value;
21042145 dv_type typ;
2146+ double relax=0.0 ;
21052147 for (int i=0 ;i<_nDecisionVars;i++)
21062148 {
21072149 p = _pDecisionVars[i]->p_index ;
21082150 value = _pDecisionVars[i]->value ;
21092151 typ = _pDecisionVars[i]->dvar_type ;
21102152
2111- if (typ == DV_STAGE) {h_iter[p]= value; }
2112- else if (typ == DV_QOUT ) {Q_iter[p]= value; }
2153+ if (typ == DV_STAGE) {ahiterHist[p][iter]= h_iter[p]; h_iter[p]=( 1.0 -relax)* value+(relax)*h_iter[p]; }
2154+ else if (typ == DV_QOUT ) {aQiterHist[p][iter]= Q_iter[p]; Q_iter[p]=( 1.0 -relax)* value+(relax)*h_iter[p]; }
21132155 }
21142156
2157+ /* for (int p = 0; p < pModel->GetNumSubBasins(); p++) {
2158+ if (pModel->GetSubBasin(p)->GetReservoir()!=NULL){
2159+ long long SBID=pModel->GetSubBasin(p)->GetID();
2160+ cout << " >Q " << SBID << ": "; for (int i = 0; i <= iter; i++) { cout << aQiterHist[p][i] << " "; }cout << endl;
2161+ //break;
2162+ }
2163+ }*/
2164+
21152165 // Calculate Non-linear solver residual
21162166 // sum of squared relative changes
21172167 // ----------------------------------------------------------------
0 commit comments