ROL
example_05.cpp
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 #include "example_05.hpp"
11 
12 typedef double RealT;
13 
14 template<class Real>
15 Real random(const ROL::Ptr<const Teuchos::Comm<int> > &comm) {
16  Real val = 0.0;
17  if ( Teuchos::rank<int>(*comm)==0 ) {
18  val = (Real)rand()/(Real)RAND_MAX;
19  }
20  Teuchos::broadcast<int,Real>(*comm,0,1,&val);
21  return val;
22 }
23 
24 int main(int argc, char* argv[]) {
25 
26  Teuchos::GlobalMPISession mpiSession(&argc, &argv);
27  ROL::Ptr<const Teuchos::Comm<int> > comm
28  = ROL::toPtr(Teuchos::DefaultComm<int>::getComm());
29 
30  // This little trick lets us print to std::cout only if a (dummy) command-line argument is provided.
31  int iprint = argc - 1;
32  ROL::Ptr<std::ostream> outStream;
33  ROL::nullstream bhs; // outputs nothing
34  if (iprint > 0 && Teuchos::rank<int>(*comm)==0)
35  outStream = ROL::makePtrFromRef(std::cout);
36  else
37  outStream = ROL::makePtrFromRef(bhs);
38 
39  int errorFlag = 0;
40 
41  try {
42  /**********************************************************************************************/
43  /************************* CONSTRUCT ROL ALGORITHM ********************************************/
44  /**********************************************************************************************/
45  // Get ROL parameterlist
46  std::string filename = "input.xml";
47  auto parlist = ROL::getParametersFromXmlFile( filename );
48  // Build ROL algorithm
49  parlist->sublist("Status Test").set("Gradient Tolerance",1.e-7);
50  parlist->sublist("Status Test").set("Step Tolerance",1.e-14);
51  parlist->sublist("Status Test").set("Iteration Limit",100);
52  /**********************************************************************************************/
53  /************************* CONSTRUCT VECTORS **************************************************/
54  /**********************************************************************************************/
55  // Build control vectors
56  int nx = 256;
57  // Construct storage for optimal solution
58  ROL::Ptr<std::vector<RealT> > z_ptr = ROL::makePtr<std::vector<RealT>>(nx+2,0);
59  ROL::Ptr<ROL::Vector<RealT> > zp = ROL::makePtr<ROL::StdVector<RealT>>(z_ptr);
60  ROL::Ptr<std::vector<RealT> > x1_ptr = ROL::makePtr<std::vector<RealT>>(nx+2,0);
61  ROL::Ptr<ROL::Vector<RealT> > x1p = ROL::makePtr<ROL::StdVector<RealT>>(x1_ptr);
62  ROL::Ptr<std::vector<RealT> > x2_ptr = ROL::makePtr<std::vector<RealT>>(nx+2,0);
63  ROL::Ptr<ROL::Vector<RealT> > x2p = ROL::makePtr<ROL::StdVector<RealT>>(x2_ptr);
64  ROL::Ptr<std::vector<RealT> > x3_ptr = ROL::makePtr<std::vector<RealT>>(nx+2,0);
65  ROL::Ptr<ROL::Vector<RealT> > x3p = ROL::makePtr<ROL::StdVector<RealT>>(x3_ptr);
66  std::vector<ROL::Ptr<ROL::Vector<RealT> > > xvec = {x1p, x2p, x3p};
67  // Create vectors for derivative check
68  ROL::Ptr<std::vector<RealT> > xr_ptr = ROL::makePtr<std::vector<RealT>>(nx+2,0);
69  ROL::StdVector<RealT> xr(xr_ptr);
70  ROL::Ptr<std::vector<RealT> > d_ptr = ROL::makePtr<std::vector<RealT>>(nx+2,0);
71  ROL::StdVector<RealT> d(d_ptr);
72  for ( int i = 0; i < nx+2; i++ ) {
73  (*xr_ptr)[i] = random<RealT>(comm);
74  (*d_ptr)[i] = random<RealT>(comm);
75  }
76  // Build state and adjoint vectors
77  ROL::Ptr<std::vector<RealT> > u_ptr = ROL::makePtr<std::vector<RealT>>(nx,1);
78  ROL::Ptr<ROL::Vector<RealT> > up = ROL::makePtr<ROL::StdVector<RealT>>(u_ptr);
79  ROL::Ptr<std::vector<RealT> > p_ptr = ROL::makePtr<std::vector<RealT>>(nx,0);
80  ROL::Ptr<ROL::Vector<RealT> > pp = ROL::makePtr<ROL::StdVector<RealT>>(p_ptr);
81  /**********************************************************************************************/
82  /************************* CONSTRUCT SOL COMPONENTS *******************************************/
83  /**********************************************************************************************/
84  // Build samplers
85  int dim = 4, nSamp = 100;
86  std::vector<RealT> tmp = {-1, 1};
87  std::vector<std::vector<RealT> > bounds(dim,tmp);
88  ROL::Ptr<ROL::BatchManager<RealT> > bman
89  = ROL::makePtr<ROL::StdTeuchosBatchManager<RealT,int>>(comm);
90  ROL::Ptr<ROL::SampleGenerator<RealT> > sampler
91  = ROL::makePtr<ROL::MonteCarloGenerator<RealT>>(nSamp,bounds,bman,false,false,100);
92  /**********************************************************************************************/
93  /************************* CONSTRUCT OBJECTIVE FUNCTION ***************************************/
94  /**********************************************************************************************/
95  // Build risk-averse objective function
96  RealT alpha = 1.e-3;
97  ROL::Ptr<ROL::Objective_SimOpt<RealT> > pobjSimOpt
98  = ROL::makePtr<Objective_BurgersControl<RealT>>(alpha,nx);
99  ROL::Ptr<ROL::Constraint_SimOpt<RealT> > pconSimOpt
100  = ROL::makePtr<Constraint_BurgersControl<RealT>>(nx);
101  pconSimOpt->setSolveParameters(*parlist);
102  ROL::Ptr<ROL::Objective<RealT> > pObj
103  = ROL::makePtr<ROL::Reduced_Objective_SimOpt<RealT>>(pobjSimOpt,pconSimOpt,up,zp,pp);
104  // Test parametrized objective functions
105  *outStream << "Check Derivatives of Parametrized Objective Function\n";
106  xvec[0]->set(xr);
107  pObj->setParameter(sampler->getMyPoint(0));
108  pObj->checkGradient(*xvec[0],d,true,*outStream);
109  pObj->checkHessVec(*xvec[0],d,true,*outStream);
110  /**********************************************************************************************/
111  /************************* SMOOTHED CVAR 1.e-2, 1.e-4, 1.e-6 **********************************/
112  /**********************************************************************************************/
113  const RealT cl(0.9), cc(1), lb(-0.5), ub(0.5);
114  const std::string ra = "Risk Averse", rm = "CVaR", dist = "Parabolic";
115  const bool storage = true;
116  RealT eps(1.e-2);
117  std::vector<RealT> stat(3,0);
118  ROL::Ptr<ROL::OptimizationProblem<RealT>> optProb;
119  ROL::Ptr<ROL::OptimizationSolver<RealT>> solver;
120  for (int i = 0; i < 3; ++i) {
121  *outStream << "\nSOLVE SMOOTHED CONDITIONAL VALUE AT RISK WITH TRUST REGION\n";
122  // Build CVaR risk measure
123  ROL::ParameterList list;
124  list.sublist("SOL").set("Type",ra);
125  list.sublist("SOL").set("Store Sampled Value and Gradient",storage);
126  list.sublist("SOL").sublist("Risk Measure").set("Name",rm);
127  list.sublist("SOL").sublist("Risk Measure").sublist(rm).set("Confidence Level",cl);
128  list.sublist("SOL").sublist("Risk Measure").sublist(rm).set("Convex Combination Parameter",cc);
129  list.sublist("SOL").sublist("Risk Measure").sublist(rm).set("Smoothing Parameter",eps);
130  list.sublist("SOL").sublist("Risk Measure").sublist(rm).sublist("Distribution").set("Name",dist);
131  list.sublist("SOL").sublist("Risk Measure").sublist(rm).sublist("Distribution").sublist(dist).set("Lower Bound",lb);
132  list.sublist("SOL").sublist("Risk Measure").sublist(rm).sublist("Distribution").sublist(dist).set("Upper Bound",ub);
133  // Build stochastic problem
134  if ( i==0 ) { xvec[i]->zero(); }
135  else { xvec[i]->set(*xvec[i-1]); }
136  optProb = ROL::makePtr<ROL::OptimizationProblem<RealT>>(pObj,xvec[i]);
137  RealT init_stat(1);
138  if ( i > 0 ) { init_stat = stat[i-1]; }
139  list.sublist("SOL").set("Initial Statistic",init_stat);
140  optProb->setStochasticObjective(list,sampler);
141  optProb->check(*outStream);
142  // Run ROL algorithm
143  parlist->sublist("Step").set("Type","Trust Region");
144  solver = ROL::makePtr<ROL::OptimizationSolver<RealT>>(*optProb,*parlist);
145  clock_t start = clock();
146  solver->solve(*outStream);
147  *outStream << "Optimization time: " << (RealT)(clock()-start)/(RealT)CLOCKS_PER_SEC << " seconds.\n";
148  // Get solution statistic
149  stat[i] = optProb->getSolutionStatistic();
150  // Update smoothing parameter
151  eps *= static_cast<RealT>(1.e-2);
152  }
153  /**********************************************************************************************/
154  /************************* NONSMOOTH PROBLEM **************************************************/
155  /**********************************************************************************************/
156  *outStream << "\nSOLVE NONSMOOTH CVAR PROBLEM WITH BUNDLE TRUST REGION\n";
157  ROL::ParameterList list;
158  list.sublist("SOL").set("Type",ra);
159  list.sublist("SOL").set("Store Sampled Value and Gradient",storage);
160  list.sublist("SOL").sublist("Risk Measure").set("Name",rm);
161  list.sublist("SOL").sublist("Risk Measure").sublist(rm).set("Confidence Level",cl);
162  list.sublist("SOL").sublist("Risk Measure").sublist(rm).set("Convex Combination Parameter",cc);
163  list.sublist("SOL").sublist("Risk Measure").sublist(rm).set("Smoothing Parameter",0.);
164  list.sublist("SOL").sublist("Risk Measure").sublist(rm).sublist("Distribution").set("Name","Dirac");
165  list.sublist("SOL").sublist("Risk Measure").sublist(rm).sublist("Distribution").sublist("Dirac").set("Location",0.);
166  // Build stochastic problem
167  zp->set(*xvec[2]);
168  optProb = ROL::makePtr<ROL::OptimizationProblem<RealT>>(pObj,zp);
169  list.sublist("SOL").set("Initial Statistic",stat[2]);
170  optProb->setStochasticObjective(list,sampler);
171  optProb->check(*outStream);
172  // Run ROL algorithm
173  parlist->sublist("Status Test").set("Iteration Limit",1000);
174  parlist->sublist("Step").sublist("Bundle").set("Epsilon Solution Tolerance",1.e-7);
175  parlist->sublist("Step").set("Type","Bundle");
176  solver = ROL::makePtr<ROL::OptimizationSolver<RealT>>(*optProb,*parlist);
177  clock_t start = clock();
178  solver->solve(*outStream);
179  *outStream << "Optimization time: " << (RealT)(clock()-start)/(RealT)CLOCKS_PER_SEC << " seconds.\n";
180  /**********************************************************************************************/
181  /************************* COMPUTE ERROR ******************************************************/
182  /**********************************************************************************************/
183  ROL::Ptr<ROL::Vector<RealT> > cErr = zp->clone();
184  RealT zstat = optProb->getSolutionStatistic();
185  *outStream << "\nSUMMARY:\n";
186  *outStream << " ---------------------------------------------\n";
187  *outStream << " True Value-At-Risk = " << zstat << "\n";
188  *outStream << " ---------------------------------------------\n";
189  RealT VARerror = std::abs(zstat-stat[0]);
190  cErr->set(*xvec[0]); cErr->axpy(-1.0,*zp);
191  RealT CTRLerror = cErr->norm();
192  RealT TOTerror1 = std::sqrt(std::pow(VARerror,2)+std::pow(CTRLerror,2));
193  *outStream << " Value-At-Risk (1.e-2) = " << stat[0] << "\n";
194  *outStream << " Value-At-Risk Error = " << VARerror << "\n";
195  *outStream << " Control Error = " << CTRLerror << "\n";
196  *outStream << " Total Error = " << TOTerror1 << "\n";
197  *outStream << " ---------------------------------------------\n";
198  VARerror = std::abs(zstat-stat[1]);
199  cErr->set(*xvec[1]); cErr->axpy(-1.0,*zp);
200  CTRLerror = cErr->norm();
201  RealT TOTerror2 = std::sqrt(std::pow(VARerror,2)+std::pow(CTRLerror,2));
202  *outStream << " Value-At-Risk (1.e-4) = " << stat[1] << "\n";
203  *outStream << " Value-At-Risk Error = " << VARerror << "\n";
204  *outStream << " Control Error = " << CTRLerror << "\n";
205  *outStream << " Total Error = " << TOTerror2 << "\n";
206  *outStream << " ---------------------------------------------\n";
207  VARerror = std::abs(zstat-stat[2]);
208  cErr->set(*xvec[2]); cErr->axpy(-1.0,*zp);
209  CTRLerror = cErr->norm();
210  RealT TOTerror3 = std::sqrt(std::pow(VARerror,2)+std::pow(CTRLerror,2));
211  *outStream << " Value-At-Risk (1.e-6) = " << stat[2] << "\n";
212  *outStream << " Value-At-Risk Error = " << VARerror << "\n";
213  *outStream << " Control Error = " << CTRLerror << "\n";
214  *outStream << " Total Error = " << TOTerror3 << "\n";
215  *outStream << " ---------------------------------------------\n\n";
216  // Comparison
217  errorFlag += ((TOTerror1 < 90.*TOTerror2) && (TOTerror2 < 90.*TOTerror3)) ? 1 : 0;
218 
219  // Output controls
220  std::ofstream control;
221  control.open("example04_control.txt");
222  for (int n = 0; n < nx+2; n++) {
223  control << std::scientific << std::setprecision(15)
224  << std::setw(25) << static_cast<RealT>(n)/static_cast<RealT>(nx+1)
225  << std::setw(25) << (*z_ptr)[n]
226  << std::endl;
227  }
228  control.close();
229 
230  }
231  catch (std::logic_error& err) {
232  *outStream << err.what() << "\n";
233  errorFlag = -1000;
234  }; // end try
235 
236  if (errorFlag != 0)
237  std::cout << "End Result: TEST FAILED\n";
238  else
239  std::cout << "End Result: TEST PASSED\n";
240 
241  return 0;
242 }
example_05.hpp
ROL::StdVector
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Definition: ROL_StdVector.hpp:27
RealT
double RealT
Definition: function/constraint/test_01.cpp:29
random
Real random(const ROL::Ptr< const Teuchos::Comm< int > > &comm)
Definition: example_05.cpp:15
ROL::details::nullstream
basic_nullstream< char, char_traits< char >> nullstream
Definition: ROL_Stream.hpp:38
main
int main(int argc, char *argv[])
Definition: function/constraint/test_01.cpp:31
dim
constexpr auto dim
Definition: vector/test_11.cpp:23
Generated by   doxygen 1.8.5