ROL
ROL_ProgressiveHedging.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_PROGRESSIVEHEDGING_H
11 #define ROL_PROGRESSIVEHEDGING_H
12 
14 #include "ROL_Solver.hpp"
15 #include "ROL_PH_Objective.hpp"
16 #include "ROL_PH_StatusTest.hpp"
17 #include "ROL_RiskVector.hpp"
20 
51 namespace ROL {
52 
53 template <class Real>
55 private:
56  const Ptr<Problem<Real>> input_;
57  const Ptr<SampleGenerator<Real>> sampler_;
58  ParameterList parlist_;
62  Real maxPen_;
63  Real update_;
64  int freq_;
65  Real ztol_;
66  int maxit_;
67  bool print_;
68 
69  bool hasStat_;
70  Ptr<PH_Objective<Real>> ph_objective_;
71  Ptr<Vector<Real>> ph_vector_;
72  Ptr<BoundConstraint<Real>> ph_bound_;
73  Ptr<Constraint<Real>> ph_constraint_;
74  Ptr<Problem<Real>> ph_problem_;
75  Ptr<Solver<Real>> ph_solver_;
76  Ptr<PH_StatusTest<Real>> ph_status_;
77  Ptr<Vector<Real>> z_psum_, z_gsum_;
78  std::vector<Ptr<Vector<Real>>> wvec_;
79 
80  void presolve(void) {
82  for (int j = 0; j < sampler_->numMySamples(); ++j) {
83  input_->getObjective()->setParameter(sampler_->getMyPoint(j));
84  if (input_->getConstraint() != nullPtr) {
85  input_->getConstraint()->setParameter(sampler_->getMyPoint(j));
86  }
87  solver.solve();
88  z_psum_->axpy(sampler_->getMyWeight(j),*input_->getPrimalOptimizationVector());
89  solver.reset();
90  }
91  // Aggregation
92  z_gsum_->zero();
93  sampler_->sumAll(*z_psum_,*z_gsum_);
94  }
95 
96 public:
97  ProgressiveHedging(const Ptr<Problem<Real>> &input,
98  const Ptr<SampleGenerator<Real>> &sampler,
99  ParameterList &parlist)
100  : input_(input), sampler_(sampler), parlist_(parlist), hasStat_(false) {
101  // Get algorithmic parameters
102  usePresolve_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Use Presolve",true);
103  useInexact_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Use Inexact Solve",true);
104  penaltyParam_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Initial Penalty Parameter",10.0);
105  maxPen_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Maximum Penalty Parameter",-1.0);
106  update_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Penalty Update Scale",10.0);
107  freq_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Penalty Update Frequency",0);
108  ztol_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Nonanticipativity Constraint Tolerance",1e-4);
109  maxit_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Iteration Limit",100);
110  print_ = parlist.sublist("SOL").sublist("Progressive Hedging").get("Print Subproblem Solve History",false);
111  maxPen_ = (maxPen_ <= static_cast<Real>(0) ? ROL_INF<Real>() : maxPen_);
112  penaltyParam_ = std::min(penaltyParam_,maxPen_);
113  // Create progressive hedging vector
114  ParameterList olist; olist.sublist("SOL") = parlist.sublist("SOL").sublist("Objective");
115  std::string type = olist.sublist("SOL").get("Type","Risk Neutral");
116  std::string prob = olist.sublist("SOL").sublist("Probability").get("Name","bPOE");
117  hasStat_ = ((type=="Risk Averse") ||
118  (type=="Deviation") ||
119  (type=="Probability" && prob=="bPOE"));
120  Ptr<ParameterList> parlistptr = makePtrFromRef<ParameterList>(olist);
121  if (hasStat_) {
122  ph_vector_ = makePtr<RiskVector<Real>>(parlistptr,
123  input_->getPrimalOptimizationVector());
124  }
125  else {
126  ph_vector_ = input_->getPrimalOptimizationVector();
127  }
128  // Create progressive hedging objective function
129  ph_objective_ = makePtr<PH_Objective<Real>>(input_->getObjective(),
130  ph_vector_,
132  olist);
133  // Create progressive hedging bound constraint
134  if (hasStat_) {
135  ph_bound_ = makePtr<RiskBoundConstraint<Real>>(parlistptr,
136  input_->getBoundConstraint());
137  }
138  else {
139  ph_bound_ = input_->getBoundConstraint();
140  }
141  // Create progressive hedging constraint
142  ph_constraint_ = nullPtr;
143  if (input_->getConstraint() != nullPtr) {
144  if (hasStat_) {
145  ph_constraint_ = makePtr<RiskLessConstraint<Real>>(input_->getConstraint());
146  }
147  else {
148  ph_constraint_ = input_->getConstraint();
149  }
150  }
151  // Build progressive hedging subproblems
152  ph_problem_ = makePtr<Problem<Real>>(ph_objective_, ph_vector_);
153  if (ph_bound_ != nullPtr) {
154  if (ph_bound_->isActivated()) {
155  ph_problem_->addBoundConstraint(ph_bound_);
156  }
157  }
158  if (ph_constraint_ != nullPtr) {
159  ph_problem_->addConstraint("PH Constraint",ph_constraint_,
160  input_->getMultiplierVector());
161  }
162  // Build progressive hedging subproblem solver
163  ph_solver_ = makePtr<Solver<Real>>(ph_problem_, parlist);
164  // Build progressive hedging status test for inexact solves
165  if (useInexact_) {
166  ph_status_ = makePtr<PH_StatusTest<Real>>(parlist,
167  *ph_vector_);
168  }
169  else {
170  ph_status_ = nullPtr;
171  }
172  // Initialize vector storage
173  z_psum_ = ph_problem_->getPrimalOptimizationVector()->clone();
174  z_gsum_ = ph_problem_->getPrimalOptimizationVector()->clone();
175  z_gsum_->set(*ph_problem_->getPrimalOptimizationVector());
176  wvec_.resize(sampler_->numMySamples());
177  for (int i = 0; i < sampler_->numMySamples(); ++i) {
178  wvec_[i] = z_psum_->clone(); wvec_[i]->zero();
179  }
180  if (usePresolve_) {
181  presolve();
182  }
183  }
184 
185  void check(std::ostream &outStream = std::cout, const int numSamples = 1) {
186  int nsamp = std::min(sampler_->numMySamples(),numSamples);
187  for (int i = 0; i < nsamp; ++i) {
188  ph_objective_->setParameter(sampler_->getMyPoint(i));
190  if (ph_constraint_ != nullPtr) {
191  ph_constraint_->setParameter(sampler_->getMyPoint(i));
192  }
193  ph_problem_->check(true,outStream);
194  }
195  }
196 
197  void run(std::ostream &outStream = std::cout) {
198  const Real zero(0);
199  std::vector<Real> vec_p(2), vec_g(2);
200  Real znorm(ROL_INF<Real>()), zdotz(0);
201  int iter(0), tspiter(0), flag = 1;
202  bool converged = true;
203  // Output
204  outStream << std::scientific << std::setprecision(6);
205  outStream << std::endl << "Progressive Hedging"
206  << std::endl << " "
207  << std::setw(8) << std::left << "iter"
208  << std::setw(15) << std::left << "||z-Ez||"
209  << std::setw(15) << std::left << "penalty"
210  << std::setw(10) << std::left << "subiter"
211  << std::setw(8) << std::left << "success"
212  << std::endl;
213  for (iter = 0; iter < maxit_; ++iter) {
214  z_psum_->zero(); vec_p[0] = zero; vec_p[1] = zero;
215  ph_problem_->getPrimalOptimizationVector()->set(*z_gsum_);
216  // Solve concurrent optimization problems
217  for (int j = 0; j < sampler_->numMySamples(); ++j) {
219  ph_problem_->getObjective()->setParameter(sampler_->getMyPoint(j));
220  if (useInexact_) {
221  ph_status_->setData(iter,z_gsum_);
222  }
223  if (ph_problem_->getConstraint() != nullPtr) {
224  ph_problem_->getConstraint()->setParameter(sampler_->getMyPoint(j));
225  }
226  if (print_) {
227  ph_solver_->solve(outStream,ph_status_,true);
228  }
229  else {
230  ph_solver_->solve(ph_status_,true);
231  }
232  wvec_[j]->axpy(penaltyParam_,*ph_problem_->getPrimalOptimizationVector());
233  vec_p[0] += sampler_->getMyWeight(j)
234  *ph_problem_->getPrimalOptimizationVector()->dot(
235  *ph_problem_->getPrimalOptimizationVector());
236  vec_p[1] += static_cast<Real>(ph_solver_->getAlgorithmState()->iter);
237  z_psum_->axpy(sampler_->getMyWeight(j),*ph_problem_->getPrimalOptimizationVector());
238  converged = (ph_solver_->getAlgorithmState()->statusFlag == EXITSTATUS_CONVERGED
239  ||ph_solver_->getAlgorithmState()->statusFlag == EXITSTATUS_USERDEFINED
240  ? converged : false);
241  ph_solver_->reset();
242  }
243  // Aggregation
244  z_gsum_->zero(); vec_g[0] = zero; vec_g[1] = zero;
245  sampler_->sumAll(*z_psum_,*z_gsum_);
246  sampler_->sumAll(&vec_p[0],&vec_g[0],2);
247  // Multiplier Update
248  for (int j = 0; j < sampler_->numMySamples(); ++j) {
249  wvec_[j]->axpy(-penaltyParam_,*z_gsum_);
250  }
251  zdotz = z_gsum_->dot(*z_gsum_);
252  znorm = std::sqrt(std::abs(vec_g[0] - zdotz));
253  tspiter += static_cast<int>(vec_g[1]);
254  // Output
255  outStream << " "
256  << std::setw(8) << std::left << iter
257  << std::setw(15) << std::left << znorm
258  << std::setw(15) << std::left << penaltyParam_
259  << std::setw(10) << std::left << static_cast<int>(vec_g[1])
260  << std::setw(8) << std::left << converged
261  << std::endl;
262  // Check termination criterion
263  if (znorm <= ztol_ && converged) {
264  flag = 0;
265  outStream << "Converged: Nonanticipativity constraint tolerance satisfied!" << std::endl;
266  break;
267  }
268  converged = true;
269  // Update penalty parameter
270  if (freq_ > 0 && iter%freq_ == 0) {
272  }
273  penaltyParam_ = std::min(penaltyParam_,maxPen_);
274  }
275  if (hasStat_) {
276  input_->getPrimalOptimizationVector()->set(*dynamicPtrCast<RiskVector<Real>>(z_gsum_)->getVector());
277  }
278  else {
279  input_->getPrimalOptimizationVector()->set(*z_gsum_);
280  }
281  // Output reason for termination
282  if (iter >= maxit_ && flag != 0) {
283  outStream << "Maximum number of iterations exceeded" << std::endl;
284  }
285  outStream << "Total number of subproblem iterations per sample: "
286  << tspiter << " / " << sampler_->numGlobalSamples()
287  << " ~ " << static_cast<int>(std::ceil(static_cast<Real>(tspiter)/static_cast<Real>(sampler_->numGlobalSamples())))
288  << std::endl;
289  }
290 
291 }; // class ProgressiveHedging
292 
293 } // namespace ROL
294 
295 #endif
Provides a simplified interface for solving a wide range of optimization problems.
Definition: ROL_Solver.hpp:30
Ptr< BoundConstraint< Real > > ph_bound_
Ptr< Problem< Real > > ph_problem_
Ptr< Vector< Real > > ph_vector_
Ptr< Solver< Real > > ph_solver_
Ptr< Vector< Real > > z_psum_
const Ptr< Problem< Real > > input_
Provides the interface to solve a stochastic program using progressive hedging.
void check(std::ostream &outStream=std::cout, const int numSamples=1)
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_
int solve(const Ptr< StatusTest< Real >> &status=nullPtr, bool combineStatus=true)
Solve optimization problem with no iteration output.
std::vector< Ptr< Vector< Real > > > wvec_
void reset()
Reset both Algorithm and Step.
Ptr< Vector< Real > > z_gsum_
void run(std::ostream &outStream=std::cout)
ProgressiveHedging(const Ptr< Problem< Real >> &input, const Ptr< SampleGenerator< Real >> &sampler, ParameterList &parlist)
Ptr< PH_Objective< Real > > ph_objective_