ROL
function/test_13.cpp
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
54
58#include "ROL_RandomVector.hpp"
59#include "ROL_StdVector.hpp"
60#include "ROL_Bounds.hpp"
61
62#include "ROL_Stream.hpp"
63#include "Teuchos_GlobalMPISession.hpp"
64
65
66
67
68// Creates f(x) = <x,Dx>/2 - <x,b> where D is a diagonal matrix
69// If D_{ii}>0 for all i, then the minimizer is the solution to
70// the linear system x_i=b_i/d_i
71template<class Real>
72ROL::Ptr<ROL::Objective<Real>>
74 const ROL::Vector<Real> &b ) {
75
76 auto op = ROL::makePtr<ROL::DiagonalOperator<Real>>(a);
77 auto vec = b.clone();
78 vec->set(b);
79 auto obj = ROL::makePtr<ROL::QuadraticObjective<Real>>(op,vec);
80 return obj;
81}
82
83typedef double RealT;
84
85int main(int argc, char *argv[]) {
86
87// typedef ROL::Vector<RealT> V;
88 typedef ROL::StdVector<RealT> SV;
89
90
91
92 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
93
94 // This little trick lets us print to std::cout only if a
95 // (dummy) command-line argument is provided.
96 int iprint = argc - 1;
97 ROL::Ptr<std::ostream> outStream;
98 ROL::nullstream bhs; // outputs nothing
99 if (iprint > 0)
100 outStream = ROL::makePtrFromRef(std::cout);
101 else
102 outStream = ROL::makePtrFromRef(bhs);
103
104 int errorFlag = 0;
105
106 try {
107
108 auto parlist = ROL::getParametersFromXmlFile("binary_constraint.xml");
109
110 // Penalty parameter
111 RealT gamma = 1.0;
112
113 RealT INF = ROL::ROL_INF<RealT>();
114 RealT NINF = ROL::ROL_NINF<RealT>();
115
116 /*----- Optimization Vector -----*/
117
118 // Initial guess
119 auto x0_ptr = ROL::makePtr<std::vector<RealT>>(4);
120 auto x0 = ROL::makePtr<SV>(x0_ptr);
122
123 auto x = x0->clone(); x->set(*x0);
124
125 /*----- Objective Function -----*/
126
127 // Diagonal quadratic objective scaling vector
128 auto d_ptr = ROL::makePtr<std::vector<RealT>>(4);
129 auto d = ROL::makePtr<SV>(d_ptr);
130
131 // Quadratic objective offset vector
132 auto b_ptr = ROL::makePtr<std::vector<RealT>>(4);
133 auto b = ROL::makePtr<SV>(b_ptr);
134
135 // Set values for objective
136 (*b_ptr)[0] = 1.0; (*b_ptr)[1] = 1.0;
137 (*b_ptr)[2] = 1.0; (*b_ptr)[3] = 1.0;
138
139 (*d_ptr)[0] = 1.0; (*d_ptr)[1] = 2.0;
140 (*d_ptr)[2] = 4.0; (*d_ptr)[3] = 8.0;
141
142 auto obj = createDiagonalQuadraticObjective( *d, *b );
143
144 // Unconstrained minimizer: x = [ 1.0, 0.5, 0.25, 0.125 ]
145
146
147 /*----- Bound Constraint -----*/
148
149 // Lower bound vector
150 auto xl_ptr = ROL::makePtr<std::vector<RealT>>(4);
151 auto xl = ROL::makePtr<SV>(xl_ptr);
152
153 // Upper bound vector
154 auto xu_ptr = ROL::makePtr<std::vector<RealT>>(4);
155 auto xu = ROL::makePtr<SV>(xu_ptr);
156
157 // Set bounds
158 (*xl_ptr)[0] = 0.0; (*xl_ptr)[1] = 0.0;
159 (*xl_ptr)[2] = NINF; (*xl_ptr)[3] = NINF;
160
161 (*xu_ptr)[0] = 1.0; (*xu_ptr)[1] = INF;
162 (*xu_ptr)[2] = 1.0; (*xu_ptr)[3] = INF;
163
164// ROL::BoundConstraint<RealT> bnd(xl,xu);
165 auto bnd = ROL::makePtr<ROL::Bounds<RealT>>(xl,xu);
166
167 /*---- Constraint and Lagrange Multiplier -----*/
168
169 auto con = ROL::makePtr<ROL::BinaryConstraint<RealT>>( bnd, gamma );
170
171 // Lagrange multiplier
172 auto l = x->dual().clone();
173 ROL::Elementwise::Fill<RealT> FillWithOnes(1.0);
174 l->applyUnary( ROL::Elementwise::Fill<RealT>(1.0) );
175
176 // Constrained minimizer set X = { [ 0, 0, 1, 0.125 ], [ 1, 0, 1, 0.125 ] }
177
178 // Create Optimization problems
179 ROL::OptimizationProblem<RealT> problem_E( obj, x, ROL::nullPtr, con, l );
180 ROL::OptimizationProblem<RealT> problem_EB( obj, x, bnd, con, l );
181
182 // Perform linear algebra and finite difference checks
183 problem_E.check( *outStream );
184
185
186 // Solve using Composite Step where the bound is not enforced and
187 // equality constraints are satisfied asymptotically
188 parlist->sublist("Step").set("Type","Composite Step");
189 ROL::OptimizationSolver<RealT> solver_cs( problem_E, *parlist );
190 solver_cs.solve( *outStream );
191 *outStream << "\n\nFinal optimization vector:";
192 x->print(*outStream);
193
194 // Reset optimization vector and Lagrange multiplier to initial values
195 x->set(*x0); l->applyUnary(FillWithOnes);
196
197
198 // Solve using Augmented Lagrangian where the bound is enforced explicitly
199 // and equality constraints are enforced through penalization
200 parlist->sublist("Step").set("Type","Augmented Lagrangian");
201 ROL::OptimizationSolver<RealT> solver_al( problem_EB, *parlist );
202 solver_al.solve( *outStream );
203 *outStream << "\n\nFinal optimization vector:";
204 x->print(*outStream);
205
206
207 // Reset optimization vector and Lagrange multiplier to initial values
208 x->set(*x0); l->applyUnary(FillWithOnes);
209
210
211 // Solve using Moreau-Yosida where the bound is enforced through penalization
212 // and equality constraints are satisfied asymptotically
213 parlist->sublist("Step").set("Type","Moreau-Yosida Penalty");
214 ROL::OptimizationSolver<RealT> solver_my( problem_EB, *parlist );
215 solver_my.solve( *outStream );
216 *outStream << "\n\nFinal optimization vector:";
217 x->print(*outStream);
218
219
220 }
221
222 catch (std::logic_error& err) {
223 *outStream << err.what() << "\n";
224 errorFlag = -1000;
225 }; // end try
226
227 if (errorFlag != 0)
228 std::cout << "End Result: TEST FAILED\n";
229 else
230 std::cout << "End Result: TEST PASSED\n";
231
232 return 0;
233}
234
Defines a no-output stream class ROL::NullStream and a function makeStreamPtr which either wraps a re...
void check(std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
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.
Provides the ROL::Vector interface for scalar values, to be used, for example, with scalar constraint...
Defines the linear algebra or vector space interface.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
int main(int argc, char *argv[])
ROL::Ptr< ROL::Objective< Real > > createDiagonalQuadraticObjective(const ROL::Vector< Real > &a, const ROL::Vector< Real > &b)
double RealT
void RandomizeVector(Vector< Real > &x, const Real &lower=0.0, const Real &upper=1.0)
Fill a ROL::Vector with uniformly-distributed random numbers in the interval [lower,...