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-12 Stanford University and the Authors. *
10 * Authors: Jack Middleton *
11 * Contributors: *
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 "InteriorPointOptimizer.h"
25
26 using std::cout;
27 using std::endl;
28
29 #ifdef _MSC_VER
30 #pragma warning(disable:4996) // don't warn about strcat, sprintf, etc.
31 #endif
32
33 namespace SimTK {
34
applicationReturnStatusToString(int status)35 static std::string applicationReturnStatusToString(int status) {
36 switch(status) {
37 case Solve_Succeeded: return "Solve succeeded";
38 case Solved_To_Acceptable_Level: return "Solved to acceptable level";
39 case Infeasible_Problem_Detected: return "Infeasible problem detected";
40 case Search_Direction_Becomes_Too_Small: return "Search direction becomes too small";
41 case Diverging_Iterates: return "Diverging iterates";
42 case User_Requested_Stop: return "User requested stop";
43 case Maximum_Iterations_Exceeded: return "Maximum iterations exceeded";
44 case Restoration_Failed: return "Restoration failed";
45 case Error_In_Step_Computation: return "Error in step computation";
46 case Not_Enough_Degrees_Of_Freedom: return "Not enough degrees of freedom";
47 case Invalid_Problem_Definition: return "Invalid problem definition";
48 case Invalid_Option: return "Invalid option";
49 case Invalid_Number_Detected: return "Invalid number detected";
50 case Unrecoverable_Exception: return "Unrecoverable exception";
51 case NonIpopt_Exception_Thrown: return "Non-Ipopt exception thrown";
52 case Insufficient_Memory: return "Insufficient memory";
53 case Internal_Error: return "Internal error";
54 default: return "Unknown Ipopt return status";
55 }
56 }
clone() const57 Optimizer::OptimizerRep* InteriorPointOptimizer::clone() const {
58 return new InteriorPointOptimizer(*this);
59 }
60
61 // Assume by the time this constructor is called, the number of parameters and constraints has been finalized
InteriorPointOptimizer(const OptimizerSystem & sys)62 InteriorPointOptimizer::InteriorPointOptimizer( const OptimizerSystem& sys )
63 : OptimizerRep( sys ) {
64
65 int n = sys.getNumParameters();
66 int m = sys.getNumConstraints();
67
68 if( n < 1 ) {
69 const char* where = " InteriorPointOptimizer Initialization";
70 const char* szName = "dimension";
71 SimTK_THROW5(SimTK::Exception::ValueOutOfRange, szName, 1, n, INT_MAX, where);
72 }
73
74 // Initialize arrays to store multipliers -- will be used for warm starts
75 mult_x_L = new Number[n];
76 mult_x_U = new Number[n];
77 mult_g = new Number[m];
78 for(int i=0;i<n;i++) mult_x_L[i] = mult_x_U[i] = 0;
79 for(int i=0;i<m;i++) mult_g[i] = 0;
80
81 g_L = new Number[m];
82 g_U = new Number[m];
83 /* set the bounds on the equality constraint functions */
84 for(int i=0;i<sys.getNumEqualityConstraints();i++){
85 g_U[i] = g_L[i] = 0.0;
86 }
87 /* set the bounds on the inequality constraint functions */
88 for(int i=sys.getNumEqualityConstraints();i<m;i++){
89 g_U[i] = SimTK::Real(POSITIVE_INF);
90 g_L[i] = 0.0;
91 }
92
93 firstOptimization = true;
94 }
95
96
optimize(Vector & results)97 SimTK::Real InteriorPointOptimizer::optimize( Vector &results ) {
98
99 int n = getOptimizerSystem().getNumParameters();
100 int m = getOptimizerSystem().getNumConstraints();
101
102 Index index_style = 0; /* C-style; start counting of rows and column indices at 0 */
103 Index nele_hess = 0;
104 Index nele_jac = n*m; /* always assume dense */
105
106 // Parameter limits
107 Number *x_L = NULL, *x_U = NULL;
108 if( getOptimizerSystem().getHasLimits() ) {
109 getOptimizerSystem().getParameterLimits( &x_L, &x_U);
110 } else {
111 x_U = new Number[n];
112 x_L = new Number[n];
113 for(int i=0;i<n;i++) {
114 x_U[i] = SimTK::Real(POSITIVE_INF);
115 x_L[i] = SimTK::Real(NEGATIVE_INF);
116 }
117 }
118
119 SimTK::Real *x = &results[0];
120
121 IpoptProblem nlp = CreateIpoptProblem(n, x_L, x_U, m, g_L, g_U, nele_jac,
122 nele_hess, index_style, objectiveFuncWrapper, constraintFuncWrapper,
123 gradientFuncWrapper, constraintJacobianWrapper, hessianWrapper);
124
125 // If you want to verify which options are getting set in the optimizer, you can create a file ipopt.opt
126 // with "print_user_options yes", and set print_level to (at least 1). It will then print the options to the screen.
127
128 // sherm 100302: you have to set all of these tolerances to get IpOpt to change
129 // its convergence criteria; see OptimalityErrorConvergenceCheck::CheckConvergence().
130 // We'll set acceptable tolerances to the same value to disable them.
131 AddIpoptNumOption(nlp, "tol", convergenceTolerance);
132 AddIpoptNumOption(nlp, "dual_inf_tol", convergenceTolerance);
133 AddIpoptNumOption(nlp, "constr_viol_tol", constraintTolerance);
134 AddIpoptNumOption(nlp, "compl_inf_tol", convergenceTolerance);
135 AddIpoptNumOption(nlp, "acceptable_tol", convergenceTolerance);
136 AddIpoptNumOption(nlp, "acceptable_dual_inf_tol", convergenceTolerance);
137 AddIpoptNumOption(nlp, "acceptable_constr_viol_tol", constraintTolerance);
138 AddIpoptNumOption(nlp, "acceptable_compl_inf_tol", convergenceTolerance);
139
140
141 AddIpoptIntOption(nlp, "max_iter", maxIterations);
142 AddIpoptStrOption(nlp, "mu_strategy", "adaptive");
143 AddIpoptStrOption(nlp, "hessian_approximation", "limited-memory"); // needs to be limited-memory unless you have explicit hessians
144 AddIpoptIntOption(nlp, "limited_memory_max_history", limitedMemoryHistory);
145 AddIpoptIntOption(nlp, "print_level", diagnosticsLevel); // default is 4
146
147 int i;
148 static const char *advancedRealOptions[] = {
149 "compl_inf_tol",
150 "dual_inf_tol",
151 "constr_viol_tol",
152 "acceptable_tol",
153 "acceptable_compl_inf_tol",
154 "acceptable_constr_viol_tol",
155 "acceptable_dual_inf_tol",
156 "diverging_iterates_tol",
157 "barrier_tol_factor",
158 "obj_scaling_factor",
159 "nlp_scaling_max_gradient",
160 "bounds_relax_factor",
161 "recalc_y_feas_tol",
162 "mu_init",
163 "mu_max_fact",
164 "mu_max",
165 "mu_min",
166 "mu_linear_decrease_factor",
167 "mu_superlinear_decrease_factor",
168 "bound_frac",
169 "bound_push",
170 "bound_mult_init_val",
171 "constr_mult_init_max",
172 "constr_mult_init_val",
173 "warm_start_bound_push",
174 "warm_start_bound_frac",
175 "warm_start_mult_bound_push",
176 "warm_start_mult_init_max",
177 "recalc_y_feas_tol",
178 "expect_infeasible_problem_ctol",
179 "soft_resto_pderror_reduction_factor",
180 "required_infeasibility_reduction",
181 "bound_mult_reset_threshold",
182 "constr_mult_reset_threshold",
183 "max_hessian_perturbation",
184 "min_hessian_perturbation",
185 "first_hessian_perturbation",
186 "perturb_inc_fact_first",
187 "perturb_inc_fact",
188 "perturb_dec_fact",
189 "jacobian_reqularization_value",
190 "derivative_test_perturbation",
191 "derivative_test_tol",
192 0};
193 Real value;
194 for(i=0;advancedRealOptions[i];i++) {
195 if(getAdvancedRealOption(advancedRealOptions[i],value))
196 AddIpoptNumOption(nlp, advancedRealOptions[i], value);
197 }
198
199 static const std::string advancedStrOptions[] = {"nlp_scaling_method",
200 "honor_original_bounds",
201 "check_derivatives_for_naninf",
202 "mu_strategy",
203 "mu_oracle",
204 "fixed_mu_oracle",
205 "alpha_for_y",
206 "recalc_y",
207 "expect_infeasible_problem",
208 "print_options_documentation",
209 "print_user_options",
210 "start_with_resto",
211 "evaluate_orig_obj_at_resto_trial",
212 "hessian_approximation",
213 "derivative_test",
214 ""};
215 std::string svalue;
216 for(i=0;!advancedStrOptions[i].empty();i++) {
217 if(getAdvancedStrOption(advancedStrOptions[i], svalue))
218 AddIpoptStrOption(nlp, advancedStrOptions[i].c_str(), svalue.c_str());
219 }
220
221 static const char* advancedIntOptions[] = {"quality_function_max_section_steps",
222 "max_soc",
223 "watchdog_shorted_iter_trigger",
224 "watchdog_trial_iter_max",
225 "max_refinement_steps",
226 "min_refinement_steps",
227 "limited_memory_max_history",
228 "limited_memory_max_skipping",
229 "derivative_test_print_all",
230 0};
231 int ivalue;
232 for(i=0;advancedIntOptions[i];i++) {
233 if(getAdvancedIntOption(advancedIntOptions[i], ivalue))
234 AddIpoptIntOption(nlp, advancedIntOptions[i], ivalue);
235 }
236
237 // Only makes sense to do a warm start if this is not the first call to optimize() (since we need
238 // reasonable starting multiplier values)
239 bool use_warm_start=false;
240 if(getAdvancedBoolOption("warm_start", use_warm_start) && use_warm_start && !firstOptimization) {
241 AddIpoptStrOption(nlp, "warm_start_init_point", "yes");
242 AddIpoptStrOption(nlp, "warm_start_entire_iterate", "yes");
243 //AddIpoptStrOption(nlp, "warm_start_same_structure", "yes"); // couldn't get this one to work
244 }
245
246 SimTK::Real obj;
247
248 int status = IpoptSolve(nlp, x, NULL, &obj, mult_g, mult_x_L, mult_x_U, (void *)this );
249
250 FreeIpoptProblem(nlp);
251
252 // Only delete these if they aren't pointing to existing parameter limits
253 if( !getOptimizerSystem().getHasLimits() ) {
254 delete [] x_U;
255 delete [] x_L;
256 }
257
258 if(status == Solved_To_Acceptable_Level) {
259 std::cout << "Ipopt: Solved to acceptable level" << std::endl;
260 } else if (status != Solve_Succeeded) {
261 if( status != NonIpopt_Exception_Thrown) {
262 char buf[1024];
263 sprintf(buf, "Ipopt: %s (status %d)",applicationReturnStatusToString(status).c_str(),status);
264 SimTK_THROW1(SimTK::Exception::OptimizerFailed, SimTK::String(buf));
265 }
266 }
267
268 firstOptimization = false;
269
270 return(obj);
271 }
272
273
274 } // namespace SimTK
275