ROL
ROL_ProgressiveHedging.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_PROGRESSIVEHEDGING_H
45#define ROL_PROGRESSIVEHEDGING_H
46
48#include "ROL_Solver.hpp"
49#include "ROL_PH_Objective.hpp"
50#include "ROL_PH_StatusTest.hpp"
51#include "ROL_RiskVector.hpp"
54
85namespace ROL {
86
87template <class Real>
89private:
90 const Ptr<Problem<Real>> input_;
91 const Ptr<SampleGenerator<Real>> sampler_;
92 ParameterList parlist_;
96 Real maxPen_;
97 Real update_;
98 int freq_;
99 Real ztol_;
101 bool print_;
102
104 Ptr<PH_Objective<Real>> ph_objective_;
105 Ptr<Vector<Real>> ph_vector_;
106 Ptr<BoundConstraint<Real>> ph_bound_;
107 Ptr<Constraint<Real>> ph_constraint_;
108 Ptr<Problem<Real>> ph_problem_;
109 Ptr<Solver<Real>> ph_solver_;
110 Ptr<PH_StatusTest<Real>> ph_status_;
111 Ptr<Vector<Real>> z_psum_, z_gsum_;
112 std::vector<Ptr<Vector<Real>>> wvec_;
113
114 void presolve(void) {
116 for (int j = 0; j < sampler_->numMySamples(); ++j) {
117 input_->getObjective()->setParameter(sampler_->getMyPoint(j));
118 if (input_->getConstraint() != nullPtr) {
119 input_->getConstraint()->setParameter(sampler_->getMyPoint(j));
120 }
121 solver.solve();
122 z_psum_->axpy(sampler_->getMyWeight(j),*input_->getPrimalOptimizationVector());
123 solver.reset();
124 }
125 // Aggregation
126 z_gsum_->zero();
127 sampler_->sumAll(*z_psum_,*z_gsum_);
128 }
129
130public:
132 const Ptr<SampleGenerator<Real>> &sampler,
133 ParameterList &parlist)
134 : input_(input), sampler_(sampler), parlist_(parlist), hasStat_(false) {
135 // Get algorithmic parameters
136 usePresolve_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Use Presolve",true);
137 useInexact_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Use Inexact Solve",true);
138 penaltyParam_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Initial Penalty Parameter",10.0);
139 maxPen_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Maximum Penalty Parameter",-1.0);
140 update_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Penalty Update Scale",10.0);
141 freq_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Penalty Update Frequency",0);
142 ztol_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Nonanticipativity Constraint Tolerance",1e-4);
143 maxit_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Iteration Limit",100);
144 print_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Print Subproblem Solve History",false);
145 maxPen_ = (maxPen_ <= static_cast<Real>(0) ? ROL_INF<Real>() : maxPen_);
147 // Create progressive hedging vector
148 ParameterList olist; olist.sublist("SOL") = parlist.sublist("SOL").sublist("Objective");
149 std::string type = olist.sublist("SOL").get("Type","Risk Neutral");
150 std::string prob = olist.sublist("SOL").sublist("Probability").get("Name","bPOE");
151 hasStat_ = ((type=="Risk Averse") ||
152 (type=="Deviation") ||
153 (type=="Probability" && prob=="bPOE"));
154 Ptr<ParameterList> parlistptr = makePtrFromRef<ParameterList>(olist);
155 if (hasStat_) {
156 ph_vector_ = makePtr<RiskVector<Real>>(parlistptr,
157 input_->getPrimalOptimizationVector());
158 }
159 else {
160 ph_vector_ = input_->getPrimalOptimizationVector();
161 }
162 // Create progressive hedging objective function
163 ph_objective_ = makePtr<PH_Objective<Real>>(input_->getObjective(),
166 olist);
167 // Create progressive hedging bound constraint
168 if (hasStat_) {
169 ph_bound_ = makePtr<RiskBoundConstraint<Real>>(parlistptr,
170 input_->getBoundConstraint());
171 }
172 else {
173 ph_bound_ = input_->getBoundConstraint();
174 }
175 // Create progressive hedging constraint
176 ph_constraint_ = nullPtr;
177 if (input_->getConstraint() != nullPtr) {
178 if (hasStat_) {
179 ph_constraint_ = makePtr<RiskLessConstraint<Real>>(input_->getConstraint());
180 }
181 else {
182 ph_constraint_ = input_->getConstraint();
183 }
184 }
185 // Build progressive hedging subproblems
186 ph_problem_ = makePtr<Problem<Real>>(ph_objective_, ph_vector_);
187 if (ph_bound_ != nullPtr) {
188 if (ph_bound_->isActivated()) {
189 ph_problem_->addBoundConstraint(ph_bound_);
190 }
191 }
192 if (ph_constraint_ != nullPtr) {
193 ph_problem_->addConstraint("PH Constraint",ph_constraint_,
194 input_->getMultiplierVector());
195 }
196 // Build progressive hedging subproblem solver
197 ph_solver_ = makePtr<Solver<Real>>(ph_problem_, parlist);
198 // Build progressive hedging status test for inexact solves
199 if (useInexact_) {
200 ph_status_ = makePtr<PH_StatusTest<Real>>(parlist,
201 *ph_vector_);
202 }
203 else {
204 ph_status_ = nullPtr;
205 }
206 // Initialize vector storage
207 z_psum_ = ph_problem_->getPrimalOptimizationVector()->clone();
208 z_gsum_ = ph_problem_->getPrimalOptimizationVector()->clone();
209 z_gsum_->set(*ph_problem_->getPrimalOptimizationVector());
210 wvec_.resize(sampler_->numMySamples());
211 for (int i = 0; i < sampler_->numMySamples(); ++i) {
212 wvec_[i] = z_psum_->clone(); wvec_[i]->zero();
213 }
214 if (usePresolve_) {
215 presolve();
216 }
217 }
218
219 void check(std::ostream &outStream = std::cout, const int numSamples = 1) {
220 int nsamp = std::min(sampler_->numMySamples(),numSamples);
221 for (int i = 0; i < nsamp; ++i) {
222 ph_objective_->setParameter(sampler_->getMyPoint(i));
224 if (ph_constraint_ != nullPtr) {
225 ph_constraint_->setParameter(sampler_->getMyPoint(i));
226 }
227 ph_problem_->check(true,outStream);
228 }
229 }
230
231 void run(std::ostream &outStream = std::cout) {
232 const Real zero(0);
233 std::vector<Real> vec_p(2), vec_g(2);
234 Real znorm(ROL_INF<Real>()), zdotz(0);
235 int iter(0), tspiter(0), flag = 1;
236 bool converged = true;
237 // Output
238 outStream << std::scientific << std::setprecision(6);
239 outStream << std::endl << "Progressive Hedging"
240 << std::endl << " "
241 << std::setw(8) << std::left << "iter"
242 << std::setw(15) << std::left << "||z-Ez||"
243 << std::setw(15) << std::left << "penalty"
244 << std::setw(10) << std::left << "subiter"
245 << std::setw(8) << std::left << "success"
246 << std::endl;
247 for (iter = 0; iter < maxit_; ++iter) {
248 z_psum_->zero(); vec_p[0] = zero; vec_p[1] = zero;
249 ph_problem_->getPrimalOptimizationVector()->set(*z_gsum_);
250 // Solve concurrent optimization problems
251 for (int j = 0; j < sampler_->numMySamples(); ++j) {
253 ph_problem_->getObjective()->setParameter(sampler_->getMyPoint(j));
254 if (useInexact_) {
255 ph_status_->setData(iter,z_gsum_);
256 }
257 if (ph_problem_->getConstraint() != nullPtr) {
258 ph_problem_->getConstraint()->setParameter(sampler_->getMyPoint(j));
259 }
260 if (print_) {
261 ph_solver_->solve(outStream,ph_status_,true);
262 }
263 else {
264 ph_solver_->solve(ph_status_,true);
265 }
266 wvec_[j]->axpy(penaltyParam_,*ph_problem_->getPrimalOptimizationVector());
267 vec_p[0] += sampler_->getMyWeight(j)
268 *ph_problem_->getPrimalOptimizationVector()->dot(
269 *ph_problem_->getPrimalOptimizationVector());
270 vec_p[1] += static_cast<Real>(ph_solver_->getAlgorithmState()->iter);
271 z_psum_->axpy(sampler_->getMyWeight(j),*ph_problem_->getPrimalOptimizationVector());
272 converged = (ph_solver_->getAlgorithmState()->statusFlag == EXITSTATUS_CONVERGED
273 ||ph_solver_->getAlgorithmState()->statusFlag == EXITSTATUS_USERDEFINED
274 ? converged : false);
275 ph_solver_->reset();
276 }
277 // Aggregation
278 z_gsum_->zero(); vec_g[0] = zero; vec_g[1] = zero;
279 sampler_->sumAll(*z_psum_,*z_gsum_);
280 sampler_->sumAll(&vec_p[0],&vec_g[0],2);
281 // Multiplier Update
282 for (int j = 0; j < sampler_->numMySamples(); ++j) {
284 }
285 zdotz = z_gsum_->dot(*z_gsum_);
286 znorm = std::sqrt(std::abs(vec_g[0] - zdotz));
287 tspiter += static_cast<int>(vec_g[1]);
288 // Output
289 outStream << " "
290 << std::setw(8) << std::left << iter
291 << std::setw(15) << std::left << znorm
292 << std::setw(15) << std::left << penaltyParam_
293 << std::setw(10) << std::left << static_cast<int>(vec_g[1])
294 << std::setw(8) << std::left << converged
295 << std::endl;
296 // Check termination criterion
297 if (znorm <= ztol_ && converged) {
298 flag = 0;
299 outStream << "Converged: Nonanticipativity constraint tolerance satisfied!" << std::endl;
300 break;
301 }
302 converged = true;
303 // Update penalty parameter
304 if (freq_ > 0 && iter%freq_ == 0) {
306 }
308 }
309 if (hasStat_) {
310 input_->getPrimalOptimizationVector()->set(*dynamicPtrCast<RiskVector<Real>>(z_gsum_)->getVector());
311 }
312 else {
313 input_->getPrimalOptimizationVector()->set(*z_gsum_);
314 }
315 // Output reason for termination
316 if (iter >= maxit_ && flag != 0) {
317 outStream << "Maximum number of iterations exceeded" << std::endl;
318 }
319 outStream << "Total number of subproblem iterations per sample: "
320 << tspiter << " / " << sampler_->numGlobalSamples()
321 << " ~ " << static_cast<int>(std::ceil(static_cast<Real>(tspiter)/static_cast<Real>(sampler_->numGlobalSamples())))
322 << std::endl;
323 }
324
325}; // class ProgressiveHedging
326
327} // namespace ROL
328
329#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Provides the interface to solve a stochastic program using progressive hedging.
ProgressiveHedging(const Ptr< Problem< Real > > &input, const Ptr< SampleGenerator< Real > > &sampler, ParameterList &parlist)
std::vector< Ptr< Vector< Real > > > wvec_
Ptr< BoundConstraint< Real > > ph_bound_
void run(std::ostream &outStream=std::cout)
Ptr< PH_StatusTest< Real > > ph_status_
const Ptr< SampleGenerator< Real > > sampler_
const Ptr< Problem< Real > > input_
Ptr< Constraint< Real > > ph_constraint_
Ptr< PH_Objective< Real > > ph_objective_
void check(std::ostream &outStream=std::cout, const int numSamples=1)
Ptr< Problem< Real > > ph_problem_
Provides a simplified interface for solving a wide range of optimization problems.
void reset()
Reset both Algorithm and Step.
int solve(const Ptr< StatusTest< Real > > &status=nullPtr, bool combineStatus=true)
Solve optimization problem with no iteration output.
virtual void set(const Vector &x)
Set where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
@ EXITSTATUS_CONVERGED
@ EXITSTATUS_USERDEFINED