1 /* -------------------------------------------------------------------------- *
2 * Simbody(tm): SimTKmath *
3 * -------------------------------------------------------------------------- *
4 * This is part of the SimTK biosimulation toolkit originating from *
5 * Simbios, the NIH National Center for Physics-Based Simulation of *
6 * Biological Structures at Stanford, funded under the NIH Roadmap for *
7 * Medical Research, grant U54 GM072970. See https://simtk.org/home/simbody. *
8 * *
9 * Portions copyright (c) 2006-13 Stanford University and the Authors. *
10 * Authors: Chris Dembia *
11 * Contributors: Jack Middleton *
12 * *
13 * Licensed under the Apache License, Version 2.0 (the "License"); you may *
14 * not use this file except in compliance with the License. You may obtain a *
15 * copy of the License at http://www.apache.org/licenses/LICENSE-2.0. *
16 * *
17 * Unless required by applicable law or agreed to in writing, software *
18 * distributed under the License is distributed on an "AS IS" BASIS, *
19 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. *
20 * See the License for the specific language governing permissions and *
21 * limitations under the License. *
22 * -------------------------------------------------------------------------- */
23
24 #include "SimTKmath.h"
25
26 #include <iostream>
27 using std::cout;
28 using std::endl;
29 using SimTK::Vector;
30 using SimTK::Matrix;
31 using SimTK::Real;
32 using SimTK::Optimizer;
33 using SimTK::OptimizerSystem;
34
35 /* This test doesn't actually run any optimizations. It simply creates an
36 * Optimizer, using the same OptimizerSystem's (not the EXACT same) used in the
37 * tests for each of the algorithms. NOTE: Since BestAvailable never selects
38 * CFSQP, we do not test for CFSQP.
39 */
40
41 /* See IpoptTest. 'BestAvailable' will select Ipopt because this sytem has
42 * constraints.
43 */
44 class IpoptSystem : public OptimizerSystem {
45 public:
46
47
objectiveFunc(const Vector & coefficients,bool new_coefficients,Real & f) const48 int objectiveFunc( const Vector &coefficients, bool new_coefficients,
49 Real& f ) const override {
50 const Real *x;
51
52 x = &coefficients[0];
53
54 f = x[0] * x[3] * (x[0] + x[1] + x[2]) + x[2];
55 return( 0 );
56 }
57
gradientFunc(const Vector & coefficients,bool new_coefficients,Vector & gradient) const58 int gradientFunc( const Vector &coefficients, bool new_coefficients,
59 Vector &gradient ) const override{
60 const Real *x;
61
62 x = &coefficients[0];
63
64 gradient[0] = x[0] * x[3] + x[3] * (x[0] + x[1] + x[2]);
65 gradient[1] = x[0] * x[3];
66 gradient[2] = x[0] * x[3] + 1;
67 gradient[3] = x[0] * (x[0] + x[1] + x[2]);
68
69 return(0);
70
71 }
72
constraintFunc(const Vector & coefficients,bool new_coefficients,Vector & constraints) const73 int constraintFunc( const Vector &coefficients, bool new_coefficients,
74 Vector &constraints ) const override{
75 const Real *x;
76
77 x = &coefficients[0];
78 constraints[0] = x[0]*x[0] + x[1]*x[1] + x[2]*x[2] + x[3]*x[3] - 40.0;
79 constraints[1] = x[0] * x[1] * x[2] * x[3] - 25.0;
80
81 return(0);
82 }
83
constraintJacobian(const Vector & coefficients,bool new_coefficients,Matrix & jac) const84 int constraintJacobian( const Vector& coefficients, bool new_coefficients,
85 Matrix& jac ) const override{
86 const Real *x;
87
88 x = &coefficients[0];
89
90 jac(0,0) = 2*x[0];
91 jac(0,1) = 2*x[1];
92 jac(0,2) = 2*x[2];
93 jac(0,3) = 2*x[3];
94 jac(1,0) = x[1]*x[2]*x[3];
95 jac(1,1) = x[0]*x[2]*x[3];
96 jac(1,2) = x[0]*x[1]*x[3];
97 jac(1,3) = x[0]*x[1]*x[2];
98
99 return(0);
100 }
101
102
IpoptSystem()103 IpoptSystem() : OptimizerSystem( 4 )
104 {
105 setNumEqualityConstraints( 1 );
106 setNumInequalityConstraints( 1 );
107 }
108
109 };
110
111
112 /* See LBFGSBTest.cpp. 'BestAvailable' will select LBFGSB because this system
113 * does NOT have constraints and does have parameter limits.
114 */
115 class LBFGSBSystem : public OptimizerSystem {
116
117 public:
118
LBFGSBSystem()119 LBFGSBSystem() : OptimizerSystem( 25 ) {}
120
objectiveFunc(const Vector & coefficients,bool new_coefficients,Real & f) const121 int objectiveFunc( const Vector &coefficients, bool new_coefficients, Real& f ) const override {
122 int i;
123
124 const Real *x = &coefficients[0];
125
126 f = .25 *(x[0]-1.0)*(x[0]-1.0);
127 for(i=1;i<getNumParameters();i++) {
128 f = f + pow(x[i]-x[i-1]*x[i-1], 2.0);
129 }
130
131 f = 4.0* f;
132 return( 0 );
133 }
134
gradientFunc(const Vector & coefficients,bool new_coefficients,Vector & gradient) const135 int gradientFunc( const Vector &coefficients, bool new_coefficients, Vector &gradient ) const override {
136 const Real *x;
137 Real t1,t2;
138 int i;
139
140 x = &coefficients[0];
141
142 t1 = x[1]-(x[0]*x[0]);
143 gradient[0] = 2.0*(x[0]-1.0)-16.0*x[0]*t1;
144 for(i=1;i<getNumParameters()-1;i++) {
145 t2=t1;
146 t1=x[i+1]-(x[i]*x[i]);
147 gradient[i]=8.0*t2-16.0*x[i]*t1;
148 }
149 gradient[getNumParameters()-1]=8.0*t1;
150
151 return(0);
152
153 }
154
155 };
156
157
158 /* See LBFGSTest.cpp. 'BestAvailable' will select LBFGS because this system
159 * does NOT have constraints and does NOT have parameter limits.
160 */
161 class LBFGSSystem : public OptimizerSystem {
162 public:
163
LBFGSSystem()164 LBFGSSystem() : OptimizerSystem( 2 ) {}
165
objectiveFunc(const Vector & coefficients,bool new_coefficients,Real & f) const166 int objectiveFunc( const Vector &coefficients, bool new_coefficients,
167 Real& f ) const override {
168
169 const Real x = coefficients[0];
170 const Real y = coefficients[1];
171
172 f = 0.5*(3*x*x+4*x*y+6*y*y) - 2*x + 8*y;
173
174 return(0);
175
176 }
177
gradientFunc(const Vector & coefficients,bool new_coefficients,Vector & gradient) const178 int gradientFunc( const Vector &coefficients, bool new_coefficients,
179 Vector &gradient ) const override {
180
181 const Real x = coefficients[0];
182 const Real y = coefficients[1];
183
184 gradient[0] = 3*x + 2*y -2;
185 gradient[1] = 2*x + 6*y +8;
186
187 return(0);
188
189 }
190 };
191
192
main()193 int main() {
194
195 int i;
196
197 // IPOPT
198 // -----
199 /* create the system to be optimized */
200 IpoptSystem sysIPOPT;
201 Vector lowerBoundsIPOPT(sysIPOPT.getNumParameters());
202 Vector upperBoundsIPOPT(sysIPOPT.getNumParameters());
203 for(i=0;i<sysIPOPT.getNumParameters();i++) {
204 lowerBoundsIPOPT[i] = 1.0;
205 upperBoundsIPOPT[i] = 5.0;
206 }
207 sysIPOPT.setParameterLimits(lowerBoundsIPOPT, upperBoundsIPOPT);
208
209
210 // LBFGSB
211 // ------
212 LBFGSBSystem sysLBFGSB;
213 Vector lowerBoundsLBFGSB(sysLBFGSB.getNumParameters());
214 Vector upperBoundsLBFGSB(sysLBFGSB.getNumParameters());
215 for(i=0;i<sysLBFGSB.getNumParameters();i=i+2) { // even numbered
216 lowerBoundsLBFGSB[i] = 1.0;
217 upperBoundsLBFGSB[i] = 100.0;
218 }
219 for(i=1;i<sysLBFGSB.getNumParameters();i=i+2) { // odd numbered
220 lowerBoundsLBFGSB[i] = -100.0;
221 upperBoundsLBFGSB[i] = 100.0;
222 }
223 sysLBFGSB.setParameterLimits( lowerBoundsLBFGSB, upperBoundsLBFGSB );
224
225
226 // LBFSGB
227 // ------
228 LBFGSSystem sysLBFGS;
229
230
231 int returnValue = 0; // assume success
232
233 try {
234
235 // Optimizer will choose the BestAvailable algorithm.
236 Optimizer optIPOPT( sysIPOPT );
237 if (optIPOPT.getAlgorithm() != SimTK::InteriorPoint)
238 { returnValue = 1; }
239
240 Optimizer optLBFGSB( sysLBFGSB );
241 if (optLBFGSB.getAlgorithm() != SimTK::LBFGSB)
242 { returnValue = 1; }
243
244 Optimizer optLBFGS( sysLBFGS );
245 if (optLBFGS.getAlgorithm() != SimTK::LBFGS)
246 { returnValue = 1; }
247
248 }
249 catch (const std::exception& e) {
250 std::cout << e.what() << std::endl;
251 returnValue = 1; // failure
252 printf("BestAvailableTest.cpp: Caught exception \n");
253 }
254
255 return( returnValue );
256 }
257