1 // Copyright (C) 2007, 2008 International Business Machines and others.
2 // All Rights Reserved.
3 // This code is published under the Eclipse Public License.
4 //
5 // $Id$
6 //
7 // Authors:  Andreas Waechter             IBM    2007-06-04
8 //               derived from IpIpoptCalculatedQuantities.cpp
9 
10 #include "IpCGPenaltyCq.hpp"
11 #include "IpCGPenaltyData.hpp"
12 #include "IpTripletHelper.hpp"
13 
14 #ifdef HAVE_CMATH
15 # include <cmath>
16 #else
17 # ifdef HAVE_MATH_H
18 #  include <math.h>
19 # else
20 #  error "don't have header file for math"
21 # endif
22 #endif
23 
24 namespace Ipopt
25 {
26 #if COIN_IPOPT_VERBOSITY > 0
27   static const Index dbg_verbosity = 0;
28 #endif
29 
CGPenaltyCq(IpoptNLP * ip_nlp,IpoptData * ip_data,IpoptCalculatedQuantities * ip_cq)30   CGPenaltyCq::CGPenaltyCq(IpoptNLP* ip_nlp,
31                            IpoptData* ip_data,
32                            IpoptCalculatedQuantities* ip_cq)
33       :
34       ip_nlp_(ip_nlp),
35       ip_data_(ip_data),
36       ip_cq_(ip_cq),
37 
38       curr_fast_direct_deriv_penalty_function_cache_(1),
39       curr_jac_cd_norm_cache_(1),
40       curr_scaled_y_Amax_cache_(1),
41       curr_added_y_nrm2_cache_(1),
42       curr_penalty_function_cache_(2),
43       trial_penalty_function_cache_(5),
44       curr_direct_deriv_penalty_function_cache_(1),
45       curr_cg_pert_fact_cache_(1),
46 
47       initialize_called_(false)
48   {
49     DBG_START_METH("CGPenaltyCq::CGPenaltyCq", dbg_verbosity);
50     DBG_ASSERT(ip_nlp_);
51     DBG_ASSERT(ip_data_);
52     DBG_ASSERT(ip_cq_);
53   }
54 
~CGPenaltyCq()55   CGPenaltyCq::~CGPenaltyCq()
56   {}
57 
RegisterOptions(const SmartPtr<RegisteredOptions> & roptions)58   void CGPenaltyCq::RegisterOptions(const SmartPtr<RegisteredOptions>& roptions)
59   {}
60 
Initialize(const Journalist & jnlst,const OptionsList & options,const std::string & prefix)61   bool CGPenaltyCq::Initialize(const Journalist& jnlst,
62                                const OptionsList& options,
63                                const std::string& prefix)
64   {
65     initialize_called_ = true;
66     return true;
67   }
68 
69   //////////////////////////////////////////////////////
70   //   Methods for the Chen-Goldfarb penalty function //
71   //////////////////////////////////////////////////////
72 
73   Number
curr_jac_cd_norm(Index nrm_type)74   CGPenaltyCq::curr_jac_cd_norm(Index nrm_type)
75   {
76     DBG_START_METH("CGPenaltyCq::curr_jac_cd_norm()",
77                    dbg_verbosity);
78 
79     Number result;
80     SmartPtr<const Matrix> jac_c = ip_cq_->curr_jac_c();
81     Index nnz = TripletHelper::GetNumberEntries(*jac_c);
82     Number* values = new Number[nnz];
83     TripletHelper::FillValues(nnz, *jac_c, values);
84     Index count = 1;
85     result = 0.;
86     for (Index i=1; i<nnz; i++) {
87       if (nrm_type == 3) {
88         result = Max(result, fabs(values[i]));
89       }
90       if (nrm_type == 1) {
91         result += fabs(values[i]);
92         count ++;
93       }
94     }
95     delete [] values;
96     SmartPtr<const Matrix> jac_d = ip_cq_->curr_jac_d();
97     nnz = TripletHelper::GetNumberEntries(*jac_d);
98     values = new Number[nnz];
99     TripletHelper::FillValues(nnz, *jac_d, values);
100     for (Index i=1; i<nnz; i++) {
101       if (nrm_type == 3) {
102         result = Max(result, fabs(values[i]));
103       }
104       if (nrm_type == 1) {
105         result += fabs(values[i]);
106         count ++;
107       }
108     }
109     delete [] values;
110     if (nrm_type == 1) {
111       result = result/count;
112     }
113     return result;
114   }
115 
116   Number
curr_penalty_function()117   CGPenaltyCq::curr_penalty_function()
118   {
119     DBG_START_METH("CGPenaltyCq::curr_penalty_function()",
120                    dbg_verbosity);
121     Number result;
122     SmartPtr<const Vector> x = ip_data_->curr()->x();
123     SmartPtr<const Vector> s = ip_data_->curr()->s();
124     std::vector<const TaggedObject*> tdeps(2);
125     tdeps[0] = GetRawPtr(x);
126     tdeps[1] = GetRawPtr(s);
127     Number mu = ip_data_->curr_mu();
128     Number penalty = CGPenData().curr_penalty();
129     std::vector<Number> sdeps(2);
130     sdeps[0] = mu;
131     sdeps[1] = penalty;
132     if (!curr_penalty_function_cache_.GetCachedResult(result, tdeps, sdeps)) {
133       if (!trial_penalty_function_cache_.GetCachedResult(result, tdeps, sdeps)) {
134         result = ip_cq_->curr_barrier_obj();
135         result += penalty*ip_cq_->curr_primal_infeasibility(NORM_2);
136       }
137       curr_penalty_function_cache_.AddCachedResult(result, tdeps, sdeps);
138     }
139     DBG_ASSERT(IsFiniteNumber(result));
140     return result;
141   }
142 
143   Number
trial_penalty_function()144   CGPenaltyCq::trial_penalty_function()
145   {
146     DBG_START_METH("CGPenaltyCq::trial_penalty_function()",
147                    dbg_verbosity);
148     Number result;
149     SmartPtr<const Vector> x = ip_data_->trial()->x();
150     SmartPtr<const Vector> s = ip_data_->trial()->s();
151     std::vector<const TaggedObject*> tdeps(2);
152     tdeps[0] = GetRawPtr(x);
153     tdeps[1] = GetRawPtr(s);
154     Number mu = ip_data_->curr_mu();
155     Number penalty = CGPenData().curr_penalty();
156     std::vector<Number> sdeps(2);
157     sdeps[0] = mu;
158     sdeps[1] = penalty;
159     if (!trial_penalty_function_cache_.GetCachedResult(result, tdeps, sdeps)) {
160       if (!curr_penalty_function_cache_.GetCachedResult(result, tdeps, sdeps)) {
161         result = ip_cq_->trial_barrier_obj();
162         result += penalty*ip_cq_->trial_primal_infeasibility(NORM_2);
163       }
164       trial_penalty_function_cache_.AddCachedResult(result, tdeps, sdeps);
165     }
166     DBG_ASSERT(IsFiniteNumber(result));
167     return result;
168   }
169 
170 
171 
curr_direct_deriv_penalty_function()172   Number CGPenaltyCq::curr_direct_deriv_penalty_function()
173   {
174     DBG_START_METH("CGPenaltyCq::curr_direct_deriv_penalty_function()",
175                    dbg_verbosity);
176 
177     Number result;
178     SmartPtr<const Vector> x = ip_data_->curr()->x();
179     SmartPtr<const Vector> s = ip_data_->curr()->s();
180     SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
181     SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
182     SmartPtr<const Vector> dy_c = CGPenData().delta_cgpen()->y_c();
183     SmartPtr<const Vector> dy_d = CGPenData().delta_cgpen()->y_d();
184     SmartPtr<const Vector> dx = CGPenData().delta_cgpen()->x();
185     SmartPtr<const Vector> ds = CGPenData().delta_cgpen()->s();
186     std::vector<const TaggedObject*> tdeps(8);
187     tdeps[0] = GetRawPtr(x);
188     tdeps[1] = GetRawPtr(s);
189     tdeps[2] = GetRawPtr(y_c);
190     tdeps[3] = GetRawPtr(y_d);
191     tdeps[4] = GetRawPtr(dy_c);
192     tdeps[5] = GetRawPtr(dy_d);
193     tdeps[6] = GetRawPtr(dx);
194     tdeps[7] = GetRawPtr(ds);
195     Number mu = ip_data_->curr_mu();
196     Number penalty = CGPenData().curr_penalty();
197     std::vector<Number> sdeps(2);
198     sdeps[0] = mu;
199     sdeps[1] = penalty;
200     if (!curr_direct_deriv_penalty_function_cache_.GetCachedResult(result, tdeps, sdeps)) {
201       result = ip_cq_->curr_grad_barrier_obj_x()->Dot(*dx) +
202                ip_cq_->curr_grad_barrier_obj_s()->Dot(*ds);
203       Number curr_inf = ip_cq_->curr_primal_infeasibility(NORM_2);
204       result -= penalty*curr_inf;
205       if (curr_inf != 0.) {
206         Number fac = penalty*CGPenData().CurrPenaltyPert()/curr_inf;
207         SmartPtr<const Vector> c = ip_cq_->curr_c();
208         SmartPtr<const Vector> d_minus_s = ip_cq_->curr_d_minus_s();
209         Number result1 = c->Dot(*y_c);
210         result1 += c->Dot(*dy_c);
211         result1 += d_minus_s->Dot(*y_d);
212         result1 += d_minus_s->Dot(*dy_d);
213         result1 *= fac;
214         result += result1;
215       }
216       curr_direct_deriv_penalty_function_cache_.AddCachedResult(result, tdeps, sdeps);
217     }
218     return result;
219   }
220 
curr_fast_direct_deriv_penalty_function()221   Number CGPenaltyCq::curr_fast_direct_deriv_penalty_function()
222   {
223     DBG_START_METH("CGPenaltyCq::curr_fast_direct_deriv_penalty_function()",
224                    dbg_verbosity);
225 
226     Number result;
227     SmartPtr<const Vector> x = ip_data_->curr()->x();
228     SmartPtr<const Vector> s = ip_data_->curr()->s();
229     DBG_ASSERT(CGPenData().HaveCgPenDeltas());
230     SmartPtr<const Vector> dy_c = CGPenData().delta_cgfast()->y_c();
231     SmartPtr<const Vector> dy_d = CGPenData().delta_cgfast()->y_d();
232     SmartPtr<const Vector> dx = CGPenData().delta_cgfast()->x();
233     SmartPtr<const Vector> ds = CGPenData().delta_cgfast()->s();
234     std::vector<const TaggedObject*> tdeps(6);
235     tdeps[0] = GetRawPtr(x);
236     tdeps[1] = GetRawPtr(s);
237     tdeps[2] = GetRawPtr(dy_c);
238     tdeps[3] = GetRawPtr(dy_d);
239     tdeps[4] = GetRawPtr(dx);
240     tdeps[5] = GetRawPtr(ds);
241     Number mu = ip_data_->curr_mu();
242     Number penalty = CGPenData().curr_penalty();
243     std::vector<Number> sdeps(2);
244     sdeps[0] = mu;
245     sdeps[1] = penalty;
246     if (!curr_fast_direct_deriv_penalty_function_cache_.GetCachedResult(result, tdeps, sdeps)) {
247       result = ip_cq_->curr_grad_barrier_obj_x()->Dot(*dx) +
248                ip_cq_->curr_grad_barrier_obj_s()->Dot(*ds);
249       Number curr_inf = ip_cq_->curr_primal_infeasibility(NORM_2);
250       result -= penalty*curr_inf;
251       if (curr_inf != 0.) {
252         Number fac = penalty*CGPenData().CurrPenaltyPert()/curr_inf;
253         SmartPtr<const Vector> c = ip_cq_->curr_c();
254         SmartPtr<const Vector> d_minus_s = ip_cq_->curr_d_minus_s();
255         Number result1 = c->Dot(*dy_c);
256         result1 += d_minus_s->Dot(*dy_d);
257         result1 *= fac;
258         result += result1;
259       }
260       curr_fast_direct_deriv_penalty_function_cache_.AddCachedResult(result, tdeps, sdeps);
261     }
262     return result;
263   }
264 
curr_cg_pert_fact()265   Number CGPenaltyCq::curr_cg_pert_fact()
266   {
267     DBG_START_METH("IpoptCalculatedQuantities::curr_cg_pert_fact()",
268                    dbg_verbosity);
269 
270     Number result;
271     SmartPtr<const Vector> x = ip_data_->curr()->x();
272     SmartPtr<const Vector> s = ip_data_->curr()->s();
273     std::vector<const TaggedObject*> tdeps(2);
274     tdeps[0] = GetRawPtr(x);
275     tdeps[1] = GetRawPtr(s);
276     Number penalty = CGPenData().curr_kkt_penalty();
277     std::vector<Number> sdeps(1);
278     sdeps[0] = penalty;
279     DBG_ASSERT(penalty>0.);
280     if (!curr_cg_pert_fact_cache_.GetCachedResult(result, tdeps, sdeps)) {
281       Number eq_2nrm = ip_cq_->curr_primal_infeasibility(NORM_2);
282       result = eq_2nrm/penalty;
283       curr_cg_pert_fact_cache_.AddCachedResult(result, tdeps, sdeps);
284     }
285     return result;
286   }
287 
dT_times_barH_times_d()288   Number CGPenaltyCq::dT_times_barH_times_d()
289   {
290     DBG_START_METH("IpoptCalculatedQuantities::dT_times_barH_times_d()",
291                    dbg_verbosity);
292     Number result;
293     SmartPtr<const Vector> d_x = CGPenData().delta_cgfast()->x();
294     SmartPtr<const Vector> d_s = CGPenData().delta_cgfast()->s();
295     SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
296     SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
297     SmartPtr<const Vector> dy_c = CGPenData().delta_cgfast()->y_c();
298     SmartPtr<const Vector> dy_d = CGPenData().delta_cgfast()->y_d();
299     SmartPtr<const Vector> c = ip_cq_->curr_c();
300     SmartPtr<const Vector> d_minus_s = ip_cq_->curr_d_minus_s();
301     Number deriv_barrier_dx = ip_cq_->curr_grad_barrier_obj_x()->Dot(*d_x);
302     Number deriv_barrier_dx_ds = deriv_barrier_dx + ip_cq_->curr_grad_barrier_obj_s()->Dot(*d_s);
303     Number penalty = CGPenData().curr_penalty();
304     result = -y_c->Dot(*dy_c);
305     result -= y_d->Dot(*dy_d);
306     result *= curr_cg_pert_fact();
307     result -= deriv_barrier_dx_ds;
308     result += c->Dot(*y_c);
309     result += d_minus_s->Dot(*y_d);
310     result -= c->Dot(*dy_c);
311     result -= d_minus_s->Dot(*dy_d);
312     result += penalty*ip_cq_->curr_primal_infeasibility(NORM_2);
313 
314     return result;
315   }
316 
317 
compute_curr_cg_penalty(const Number pen_des_fact)318   Number CGPenaltyCq::compute_curr_cg_penalty(const Number pen_des_fact )
319   {
320     DBG_START_METH("CGPenaltyCq::compute_curr_cg_penalty()",
321                    dbg_verbosity);
322 
323     SmartPtr<const Vector> d_x = ip_data_->delta()->x();
324     SmartPtr<const Vector> d_s = ip_data_->delta()->s();
325     SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
326     SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
327     SmartPtr<const Vector> dy_c = ip_data_->delta()->y_c();
328     SmartPtr<const Vector> dy_d = ip_data_->delta()->y_d();
329 
330     // Compute delta barrier times (delta x, delta s)
331     Number deriv_barrier_dx = ip_cq_->curr_grad_barrier_obj_x()->Dot(*d_x);
332     Number deriv_barrier_dx_ds = deriv_barrier_dx + ip_cq_->curr_grad_barrier_obj_s()->Dot(*d_s);
333     // Compute delta x times the damped Hessian times delta x
334     SmartPtr<const Vector> tem_jac_cT_times_y_c =
335       ip_cq_->curr_jac_cT_times_vec(*y_c);
336     SmartPtr<const Vector> tem_jac_cT_times_dy_c =
337       ip_cq_->curr_jac_cT_times_vec(*dy_c);
338     SmartPtr<Vector> tem_jac_cT_times_y_c_plus_dy_c =
339       tem_jac_cT_times_y_c->MakeNew();
340     tem_jac_cT_times_y_c_plus_dy_c->AddTwoVectors(1.,*tem_jac_cT_times_y_c, 1.,
341         *tem_jac_cT_times_dy_c, 0.);
342     SmartPtr<const Vector> tem_jac_dT_times_y_d =
343       ip_cq_->curr_jac_dT_times_vec(*y_d);
344     SmartPtr<const Vector> tem_jac_dT_times_dy_d =
345       ip_cq_->curr_jac_cT_times_vec(*dy_c);
346     SmartPtr<Vector> tem_jac_dT_times_y_d_plus_dy_d =
347       tem_jac_cT_times_y_c->MakeNew();
348     tem_jac_dT_times_y_d_plus_dy_d->AddTwoVectors(1.,*tem_jac_dT_times_y_d, 1.,
349         *tem_jac_dT_times_dy_d, 0.);
350     Number d_xs_times_damped_Hessian_times_d_xs = -deriv_barrier_dx_ds;
351     d_xs_times_damped_Hessian_times_d_xs +=
352       -(tem_jac_cT_times_y_c_plus_dy_c->Dot(*d_x)
353         +tem_jac_dT_times_y_d_plus_dy_d->Dot(*d_x)
354         -y_d->Dot(*d_s)
355         -dy_d->Dot(*d_s));
356     Number dxs_nrm = pow(d_x->Nrm2(), 2.) + pow(d_s->Nrm2(), 2.);
357     d_xs_times_damped_Hessian_times_d_xs = Max(1e-8*dxs_nrm,
358                                            d_xs_times_damped_Hessian_times_d_xs);
359     Number infeasibility = ip_cq_->curr_primal_infeasibility(NORM_2);
360     Number penalty = 0.;
361     if (infeasibility > 0.) {
362       Number deriv_inf = 0.;
363       Number fac = CGPenData().CurrPenaltyPert()/infeasibility;
364       SmartPtr<const Vector> c = ip_cq_->curr_c();
365       SmartPtr<const Vector> d_minus_s = ip_cq_->curr_d_minus_s();
366       if (CGPenData().HaveCgFastDeltas()) {
367         SmartPtr<const Vector> fast_dy_c = CGPenData().delta_cgfast()->y_c();
368         SmartPtr<const Vector> fast_dy_d = CGPenData().delta_cgfast()->y_d();
369         deriv_inf += c->Dot(*fast_dy_c);
370         deriv_inf += d_minus_s->Dot(*fast_dy_d);
371         deriv_inf *= fac;
372         deriv_inf -= infeasibility;
373       }
374       else {
375         SmartPtr<const Vector> cgpen_dy_c = CGPenData().delta_cgpen()->y_c();
376         SmartPtr<const Vector> cgpen_dy_d = CGPenData().delta_cgpen()->y_d();
377         deriv_inf += c->Dot(*cgpen_dy_c);
378         deriv_inf += c->Dot(*y_c);
379         deriv_inf += d_minus_s->Dot(*cgpen_dy_d);
380         deriv_inf += d_minus_s->Dot(*y_d);
381         deriv_inf *= fac;
382         deriv_inf -= infeasibility;
383       }
384       penalty = -(deriv_barrier_dx_ds + pen_des_fact*
385                   d_xs_times_damped_Hessian_times_d_xs)/
386                 (deriv_inf + pen_des_fact*infeasibility);
387     }
388 
389     return penalty;
390   }
391 
392 
compute_curr_cg_penalty_scale()393   Number CGPenaltyCq::compute_curr_cg_penalty_scale()
394   {
395     DBG_START_METH("CGPenaltyCq::compute_curr_cg_penalty_scale()",
396                    dbg_verbosity);
397     Number penalty;
398     Number infeasibility = ip_cq_->curr_primal_infeasibility(NORM_2);
399     if (!CGPenData().NeverTryPureNewton()) {
400       penalty = Min(1e13, infeasibility*1e9);
401     }
402     else {
403       Number reference = (curr_jac_cd_norm(1) +
404                           ip_cq_->curr_primal_infeasibility(NORM_1)/
405                           (ip_data_->curr()->y_c()->Dim()+
406                            ip_data_->curr()->y_d()->Dim()))/2.;
407       if (CGPenData().restor_iter() == ip_data_->iter_count() ||
408           ip_data_->iter_count() == 0) {
409         reference_infeasibility_ = Min(1., infeasibility);
410       }
411       Number i = CGPenData().restor_counter();
412       Number fac = 4*1e-2*pow(1e1, i);
413       //Number fac = 1e-2;
414       penalty = Min(1e4,infeasibility) / (reference*fac*
415                                           pow(reference_infeasibility_, 1));
416     }
417 
418     return penalty;
419   }
420 
curr_scaled_y_Amax()421   Number CGPenaltyCq::curr_scaled_y_Amax()
422   {
423     DBG_START_METH("CGPenaltyCq::curr_scaled_y_Amax()",
424                    dbg_verbosity);
425 
426     Number result;
427     SmartPtr<const Vector> x = ip_data_->curr()->x();
428     SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
429     SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
430     std::vector<const TaggedObject*> deps(3);
431     deps[0] = GetRawPtr(x);
432     deps[1] = GetRawPtr(y_c);
433     deps[2] = GetRawPtr(y_d);
434     if (!curr_scaled_y_Amax_cache_.GetCachedResult(result, deps)) {
435       result = Max(y_c->Amax(), y_d->Amax());
436       result /= Max(1., ip_cq_->curr_grad_f()->Amax());
437       curr_scaled_y_Amax_cache_.AddCachedResult(result, deps);
438     }
439     return result;
440   }
441 
curr_added_y_nrm2()442   Number CGPenaltyCq::curr_added_y_nrm2()
443   {
444     DBG_START_METH("CGPenaltyCq::curr_added_y_nrm2()",
445                    dbg_verbosity);
446 
447     Number result;
448     SmartPtr<const Vector> x = ip_data_->curr()->x();
449     SmartPtr<const Vector> y_c = ip_data_->curr()->y_c();
450     SmartPtr<const Vector> y_d = ip_data_->curr()->y_d();
451     std::vector<const TaggedObject*> deps(3);
452     deps[0] = GetRawPtr(x);
453     deps[1] = GetRawPtr(y_c);
454     deps[2] = GetRawPtr(y_d);
455     if (!curr_added_y_nrm2_cache_.GetCachedResult(result, deps)) {
456       SmartPtr<Vector> y_c_plus_dy_c = ip_data_->delta()->y_c()->MakeNew();
457       SmartPtr<Vector> y_d_plus_dy_d = ip_data_->delta()->y_d()->MakeNew();
458       y_c_plus_dy_c->AddTwoVectors(1.,*ip_data_->delta()->y_c(),
459                                    1.,*ip_data_->curr()->y_c(), 0.);
460       y_d_plus_dy_d->AddTwoVectors(1.,*ip_data_->delta()->y_d(),
461                                    1.,*ip_data_->curr()->y_d(), 0.);
462       result = sqrt(pow(y_c_plus_dy_c->Nrm2(),2)
463                     + pow(y_d_plus_dy_d->Nrm2(),2) );
464       curr_added_y_nrm2_cache_.AddCachedResult(result, deps);
465     }
466     return result;
467   }
468 
469 } // namespace Ipopt
470