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 
85 namespace ROL {
86 
87 template <class Real>
89 private:
90  const Ptr<OptimizationProblem<Real>> input_;
91  const Ptr<SampleGenerator<Real>> sampler_;
92  ParameterList parlist_;
96  Real maxPen_;
97  Real update_;
98  int freq_;
99  Real ztol_;
100  int maxit_;
101  bool print_;
102 
103  bool hasStat_;
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<OptimizationProblem<Real>> ph_problem_;
109  Ptr<Problem<Real>> ph_problem_new_;
110  Ptr<Solver<Real>> ph_solver_;
111  Ptr<PH_StatusTest<Real>> ph_status_;
112  Ptr<Vector<Real>> z_psum_, z_gsum_;
113  std::vector<Ptr<Vector<Real>>> wvec_;
114 
115  void presolve(void) {
117  for (int j = 0; j < sampler_->numMySamples(); ++j) {
118  input_->getObjective()->setParameter(sampler_->getMyPoint(j));
119  if (input_->getConstraint() != nullPtr) {
120  input_->getConstraint()->setParameter(sampler_->getMyPoint(j));
121  }
122  solver.solve();
123  z_psum_->axpy(sampler_->getMyWeight(j),*input_->getSolutionVector());
124  solver.reset();
125  }
126  // Aggregation
127  z_gsum_->zero();
128  sampler_->sumAll(*z_psum_,*z_gsum_);
129  }
130 
131 public:
133  const Ptr<SampleGenerator<Real>> &sampler,
134  ParameterList &parlist)
135  : input_(input), sampler_(sampler), parlist_(parlist), hasStat_(false) {
136  // Get algorithmic parameters
137  usePresolve_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Use Presolve",true);
138  useInexact_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Use Inexact Solve",true);
139  penaltyParam_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Initial Penalty Parameter",10.0);
140  maxPen_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Maximum Penalty Parameter",-1.0);
141  update_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Penalty Update Scale",10.0);
142  freq_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Penalty Update Frequency",0);
143  ztol_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Nonanticipativity Constraint Tolerance",1e-4);
144  maxit_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Iteration Limit",100);
145  print_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Print Subproblem Solve History",false);
146  maxPen_ = (maxPen_ <= static_cast<Real>(0) ? ROL_INF<Real>() : maxPen_);
147  penaltyParam_ = std::min(penaltyParam_,maxPen_);
148  // Create progressive hedging vector
149  ParameterList olist; olist.sublist("SOL") = parlist.sublist("SOL").sublist("Objective");
150  std::string type = olist.sublist("SOL").get("Type","Risk Neutral");
151  std::string prob = olist.sublist("SOL").sublist("Probability").get("Name","bPOE");
152  hasStat_ = ((type=="Risk Averse") ||
153  (type=="Deviation") ||
154  (type=="Probability" && prob=="bPOE"));
155  Ptr<ParameterList> parlistptr = makePtrFromRef<ParameterList>(olist);
156  if (hasStat_) {
157  ph_vector_ = makePtr<RiskVector<Real>>(parlistptr,
158  input_->getSolutionVector());
159  }
160  else {
161  ph_vector_ = input_->getSolutionVector();
162  }
163  // Create progressive hedging objective function
164  ph_objective_ = makePtr<PH_Objective<Real>>(input_->getObjective(),
165  ph_vector_,
167  olist);
168  // Create progressive hedging bound constraint
169  if (hasStat_) {
170  ph_bound_ = makePtr<RiskBoundConstraint<Real>>(parlistptr,
171  input_->getBoundConstraint());
172  }
173  else {
174  ph_bound_ = input_->getBoundConstraint();
175  }
176  // Create progressive hedging constraint
177  ph_constraint_ = nullPtr;
178  if (input_->getConstraint() != nullPtr) {
179  if (hasStat_) {
180  ph_constraint_ = makePtr<RiskLessConstraint<Real>>(input_->getConstraint());
181  }
182  else {
183  ph_constraint_ = input_->getConstraint();
184  }
185  }
186  // Build progressive hedging subproblems
187  ph_problem_ = makePtr<OptimizationProblem<Real>>(ph_objective_,
188  ph_vector_,
189  ph_bound_,
191  input_->getMultiplierVector());
192  ph_problem_new_ = makePtr<Problem<Real>>(ph_problem_->getObjective(),
193  ph_problem_->getSolutionVector());
194  if (ph_problem_->getBoundConstraint() != nullPtr) {
195  if (ph_problem_->getBoundConstraint()->isActivated()) {
196  ph_problem_new_->addBoundConstraint(ph_problem_->getBoundConstraint());
197  }
198  }
199  if (ph_problem_->getConstraint() != nullPtr) {
200  ph_problem_new_->addConstraint("PH Constraint",ph_problem_->getConstraint(),
201  ph_problem_->getMultiplierVector());
202  }
203  // Build progressive hedging subproblem solver
204  ph_solver_ = makePtr<Solver<Real>>(ph_problem_new_, parlist);
205  // Build progressive hedging status test for inexact solves
206  if (useInexact_) {
207  ph_status_ = makePtr<PH_StatusTest<Real>>(parlist,
208  *ph_vector_);
209 //*input_->getSolutionVector());
210  }
211  else {
212  ph_status_ = nullPtr;
213  }
214  // Initialize vector storage
215  z_psum_ = ph_problem_->getSolutionVector()->clone();
216  z_gsum_ = ph_problem_->getSolutionVector()->clone();
217  z_gsum_->set(*ph_problem_->getSolutionVector());
218  wvec_.resize(sampler_->numMySamples());
219  for (int i = 0; i < sampler_->numMySamples(); ++i) {
220  wvec_[i] = z_psum_->clone(); wvec_[i]->zero();
221  }
222  if (usePresolve_) {
223  presolve();
224  }
225  }
226 
227  void check(std::ostream &outStream = std::cout, const int numSamples = 1) {
228  int nsamp = std::min(sampler_->numMySamples(),numSamples);
229  for (int i = 0; i < nsamp; ++i) {
230  ph_objective_->setParameter(sampler_->getMyPoint(i));
232  if (ph_constraint_ != nullPtr) {
233  ph_constraint_->setParameter(sampler_->getMyPoint(i));
234  }
235  ph_problem_->check(outStream);
236  }
237  }
238 
239  void run(std::ostream &outStream = std::cout) {
240  const Real zero(0);
241  std::vector<Real> vec_p(2), vec_g(2);
242  Real znorm(ROL_INF<Real>()), zdotz(0);
243  int iter(0), tspiter(0), flag = 1;
244  bool converged = true;
245  // Output
246  outStream << std::scientific << std::setprecision(6);
247  outStream << std::endl << "Progressive Hedging"
248  << std::endl << " "
249  << std::setw(8) << std::left << "iter"
250  << std::setw(15) << std::left << "||z-Ez||"
251  << std::setw(15) << std::left << "penalty"
252  << std::setw(10) << std::left << "subiter"
253  << std::setw(8) << std::left << "success"
254  << std::endl;
255  for (iter = 0; iter < maxit_; ++iter) {
256  z_psum_->zero(); vec_p[0] = zero; vec_p[1] = zero;
257  ph_problem_->getSolutionVector()->set(*z_gsum_);
258  // Solve concurrent optimization problems
259  for (int j = 0; j < sampler_->numMySamples(); ++j) {
261  ph_problem_->getObjective()->setParameter(sampler_->getMyPoint(j));
262  if (useInexact_) {
263  ph_status_->setData(iter,z_gsum_);
264  }
265  if (ph_problem_->getConstraint() != nullPtr) {
266  ph_problem_->getConstraint()->setParameter(sampler_->getMyPoint(j));
267  }
268  if (print_) {
269  ph_solver_->solve(outStream,ph_status_,true);
270  }
271  else {
272  ph_solver_->solve(ph_status_,true);
273  }
274  wvec_[j]->axpy(penaltyParam_,*ph_problem_->getSolutionVector());
275  vec_p[0] += sampler_->getMyWeight(j)
276  *ph_problem_->getSolutionVector()->dot(
277  *ph_problem_->getSolutionVector());
278  vec_p[1] += static_cast<Real>(ph_solver_->getAlgorithmState()->iter);
279  z_psum_->axpy(sampler_->getMyWeight(j),*ph_problem_->getSolutionVector());
280  converged = (ph_solver_->getAlgorithmState()->statusFlag == EXITSTATUS_CONVERGED
281  ||ph_solver_->getAlgorithmState()->statusFlag == EXITSTATUS_USERDEFINED
282  ? converged : false);
283  ph_solver_->reset();
284  }
285  // Aggregation
286  z_gsum_->zero(); vec_g[0] = zero; vec_g[1] = zero;
287  sampler_->sumAll(*z_psum_,*z_gsum_);
288  sampler_->sumAll(&vec_p[0],&vec_g[0],2);
289  // Multiplier Update
290  for (int j = 0; j < sampler_->numMySamples(); ++j) {
291  wvec_[j]->axpy(-penaltyParam_,*z_gsum_);
292  }
293  zdotz = z_gsum_->dot(*z_gsum_);
294  znorm = std::sqrt(std::abs(vec_g[0] - zdotz));
295  tspiter += static_cast<int>(vec_g[1]);
296  // Output
297  outStream << " "
298  << std::setw(8) << std::left << iter
299  << std::setw(15) << std::left << znorm
300  << std::setw(15) << std::left << penaltyParam_
301  << std::setw(10) << std::left << static_cast<int>(vec_g[1])
302  << std::setw(8) << std::left << converged
303  << std::endl;
304  // Check termination criterion
305  if (znorm <= ztol_ && converged) {
306  flag = 0;
307  outStream << "Converged: Nonanticipativity constraint tolerance satisfied!" << std::endl;
308  break;
309  }
310  converged = true;
311  // Update penalty parameter
312  if (freq_ > 0 && iter%freq_ == 0) {
314  }
315  penaltyParam_ = std::min(penaltyParam_,maxPen_);
316  }
317  if (hasStat_) {
318  input_->getSolutionVector()->set(*dynamicPtrCast<RiskVector<Real>>(z_gsum_)->getVector());
319  }
320  else {
321  input_->getSolutionVector()->set(*z_gsum_);
322  }
323  // Output reason for termination
324  if (iter >= maxit_ && flag != 0) {
325  outStream << "Maximum number of iterations exceeded" << std::endl;
326  }
327  outStream << "Total number of subproblem iterations per sample: "
328  << tspiter << " / " << sampler_->numGlobalSamples()
329  << " ~ " << static_cast<int>(std::ceil(static_cast<Real>(tspiter)/static_cast<Real>(sampler_->numGlobalSamples())))
330  << std::endl;
331  }
332 
333 }; // class ProgressiveHedging
334 
335 } // namespace ROL
336 
337 #endif
ProgressiveHedging(const Ptr< OptimizationProblem< Real >> &input, const Ptr< SampleGenerator< Real >> &sampler, ParameterList &parlist)
Ptr< BoundConstraint< Real > > ph_bound_
Ptr< Vector< Real > > ph_vector_
Ptr< Solver< Real > > ph_solver_
Provides the interface to solve a stochastic program using progressive hedging.
void check(std::ostream &outStream=std::cout, const int numSamples=1)
Ptr< OptimizationProblem< Real > > ph_problem_
Ptr< Constraint< Real > > ph_constraint_
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Ptr< PH_StatusTest< Real > > ph_status_
const Ptr< SampleGenerator< Real > > sampler_
void reset(const bool resetAlgo=true)
Reset both Algorithm and Step.
Ptr< Problem< Real > > ph_problem_new_
std::vector< Ptr< Vector< Real > > > wvec_
Provides a simplified interface for solving a wide range of optimization problems.
int solve(const ROL::Ptr< StatusTest< Real > > &status=ROL::nullPtr, const bool combineStatus=true)
Solve optimization problem with no iteration output.
const Ptr< OptimizationProblem< Real > > input_
void run(std::ostream &outStream=std::cout)
Ptr< PH_Objective< Real > > ph_objective_