1 // (C) Copyright International Business Machines Corporation 2006
2 // All Rights Reserved.
3 // This code is published under the Eclipse Public License.
4 //
5 // $Id$
6 //
7 // Authors:  Andreas Waechter          IBM    2006-03-09
8 
9 #include "BonIpoptInteriorWarmStarter.hpp"
10 #include "IpDenseVector.hpp"
11 #include "IpIpoptData.hpp"
12 #include "IpIpoptCalculatedQuantities.hpp"
13 
14 #include <cassert>
15 
16 using namespace Ipopt;
17 
18 namespace Bonmin
19 {
20   IpoptInteriorWarmStarter::
IpoptInteriorWarmStarter(Index n,const Number * x_l,const Number * x_u,Number nlp_lower_bound_inf,Number nlp_upper_bound_inf,bool store_several_iterates)21   IpoptInteriorWarmStarter(Index n,
22       const Number* x_l, const Number* x_u,
23       Number nlp_lower_bound_inf,
24       Number nlp_upper_bound_inf,
25       bool store_several_iterates)
26       :
27       nlp_lower_bound_inf_(nlp_lower_bound_inf),
28       nlp_upper_bound_inf_(nlp_upper_bound_inf),
29       store_several_iterates_(store_several_iterates),
30       n_(n),
31       n_stored_iterates_(0)
32   {
33     x_l_prev_ = new double[n];
34     x_u_prev_ = new double[n];
35 
36     for (Index i=0; i<n; i++) {
37       x_l_prev_[i] = x_l[i];
38       x_u_prev_[i] = x_u[i];
39     }
40   }
41 
42   IpoptInteriorWarmStarter::
~IpoptInteriorWarmStarter()43   ~IpoptInteriorWarmStarter()
44   {
45     delete [] x_l_prev_;
46     delete [] x_u_prev_;
47   }
48 
49   bool IpoptInteriorWarmStarter::
UpdateStoredIterates(AlgorithmMode mode,const IpoptData & ip_data,IpoptCalculatedQuantities & ip_cq)50   UpdateStoredIterates(AlgorithmMode mode,
51       const IpoptData& ip_data,
52       IpoptCalculatedQuantities& ip_cq)
53   {
54     // Don't store anything during the restoration phase
55     if (mode==RestorationPhaseMode) {
56       return true;
57     }
58 
59     // Get some useful information out of the Ipopt objects
60     Index iter = ip_data.iter_count();
61     Number mu = ip_data.curr_mu();
62     Number nlp_error = ip_cq.curr_nlp_error();
63     Number primal_inf = ip_cq.curr_primal_infeasibility(NORM_MAX);
64     Number dual_inf = ip_cq.curr_dual_infeasibility(NORM_MAX);
65     Number complementarity = ip_cq.curr_complementarity(0., NORM_MAX);
66     if (store_several_iterates_ || n_stored_iterates_==0) {
67       // For now, we just store everything
68       n_stored_iterates_++;
69       stored_iter_.push_back(iter);
70       stored_iterates_.push_back(ip_data.curr());
71       stored_mu_.push_back(mu);
72       stored_nlp_error_.push_back(nlp_error);
73       stored_primal_inf_.push_back(primal_inf);
74       stored_dual_inf_.push_back(dual_inf);
75       stored_compl_.push_back(complementarity);
76     }
77     else {
78       stored_iter_[0] = iter;
79       stored_iterates_[0] = ip_data.curr();
80       stored_mu_[0] = mu;
81       stored_nlp_error_[0] = nlp_error;
82       stored_primal_inf_[0] = primal_inf;
83       stored_dual_inf_[0] = dual_inf;
84       stored_compl_[0] = complementarity;
85     }
86     return true;
87   }
88 
89   bool IpoptInteriorWarmStarter::
Finalize()90   Finalize()
91   {
92     // For now, nothing.  Later we could clean up, reduce storage etc.
93     return true;
94   }
95 
96   bool IpoptInteriorWarmStarter::
WarmStartIterate(Index n,const Number * x_l_new,const Number * x_u_new,IteratesVector & warm_start_iterate)97   WarmStartIterate(Index n, const Number* x_l_new,
98       const Number* x_u_new,
99       IteratesVector& warm_start_iterate)
100   {
101     assert(n==n_);
102 
103     if (n_stored_iterates_==0) {
104       return false;
105     }
106 
107     // For now let's just assume that we want to restore the 4-th latest
108     // iterate from the previous solve
109 
110     Index iter_wanted = Max(0, n_stored_iterates_-5);
111 
112     SmartPtr<const Vector> prev_x = stored_iterates_[iter_wanted]->x();
113     SmartPtr<const Vector> prev_s = stored_iterates_[iter_wanted]->s();
114     SmartPtr<const Vector> prev_z_L = stored_iterates_[iter_wanted]->z_L();
115     SmartPtr<const Vector> prev_z_U = stored_iterates_[iter_wanted]->z_U();
116     SmartPtr<const Vector> prev_y_c = stored_iterates_[iter_wanted]->y_c();
117     SmartPtr<const Vector> prev_y_d = stored_iterates_[iter_wanted]->y_d();
118     SmartPtr<const Vector> prev_v_L = stored_iterates_[iter_wanted]->v_L();
119     SmartPtr<const Vector> prev_v_U = stored_iterates_[iter_wanted]->v_U();
120 
121     const DenseVector* d_x = dynamic_cast<const DenseVector*> (GetRawPtr(prev_x));
122     const DenseVector* d_s = dynamic_cast<const DenseVector*> (GetRawPtr(prev_s));
123     const DenseVector* d_z_L = dynamic_cast<const DenseVector*> (GetRawPtr(prev_z_L));
124     const DenseVector* d_z_U = dynamic_cast<const DenseVector*> (GetRawPtr(prev_z_U));
125     const DenseVector* d_y_c = dynamic_cast<const DenseVector*> (GetRawPtr(prev_y_c));
126     const DenseVector* d_y_d = dynamic_cast<const DenseVector*> (GetRawPtr(prev_y_d));
127     const DenseVector* d_v_L = dynamic_cast<const DenseVector*> (GetRawPtr(prev_v_L));
128     const DenseVector* d_v_U = dynamic_cast<const DenseVector*> (GetRawPtr(prev_v_U));
129 
130     const Number* x_vals_prev = d_x->Values();
131     const Number* s_vals_prev = d_s->Values();
132     const Number* z_L_vals_prev = d_z_L->Values();
133     const Number* z_U_vals_prev = d_z_U->Values();
134     const Number* y_c_vals_prev = d_y_c->Values();
135     const Number* y_d_vals_prev = d_y_d->Values();
136     const Number* v_L_vals_prev = d_v_L->Values();
137     const Number* v_U_vals_prev = d_v_U->Values();
138 
139     DenseVector* d_x_new = dynamic_cast<DenseVector*> (GetRawPtr(warm_start_iterate.x_NonConst()));
140     DenseVector* d_s_new = dynamic_cast<DenseVector*> (GetRawPtr(warm_start_iterate.s_NonConst()));
141     DenseVector* d_z_L_new = dynamic_cast<DenseVector*> (GetRawPtr(warm_start_iterate.z_L_NonConst()));
142     DenseVector* d_z_U_new = dynamic_cast<DenseVector*> (GetRawPtr(warm_start_iterate.z_U_NonConst()));
143     DenseVector* d_y_c_new = dynamic_cast<DenseVector*> (GetRawPtr(warm_start_iterate.y_c_NonConst()));
144     DenseVector* d_y_d_new = dynamic_cast<DenseVector*> (GetRawPtr(warm_start_iterate.y_d_NonConst()));
145     DenseVector* d_v_L_new = dynamic_cast<DenseVector*> (GetRawPtr(warm_start_iterate.v_L_NonConst()));
146     DenseVector* d_v_U_new = dynamic_cast<DenseVector*> (GetRawPtr(warm_start_iterate.v_U_NonConst()));
147 
148     Number* x_vals_new = d_x_new->Values();
149     Number* s_vals_new = d_s_new->Values();
150     Number* z_L_vals_new = d_z_L_new->Values();
151     Number* z_U_vals_new = d_z_U_new->Values();
152     Number* y_c_vals_new = d_y_c_new->Values();
153     Number* y_d_vals_new = d_y_d_new->Values();
154     Number* v_L_vals_new = d_v_L_new->Values();
155     Number* v_U_vals_new = d_v_U_new->Values();
156 
157     // Now copy the primal variables from the old to the new problem;
158     // make sure that we take care of the fixed variables
159     Index ix_prev = 0;
160     Index ix_new = 0;
161     Index izL_prev = 0;
162     Index izL_new = 0;
163     Index izU_prev = 0;
164     Index izU_new = 0;
165     for (Index i=0; i<n_; i++) {
166       if (x_l_new[i]<x_u_new[i]) {
167         DBG_ASSERT(x_l_prev_[i]<x_u_prev_[i]);
168         x_vals_new[ix_new] = x_vals_prev[ix_prev];
169         ix_new++;
170         ix_prev++;
171         if (x_l_new[i]>nlp_lower_bound_inf_) {
172           DBG_ASSERT(x_l_prev_[i]>nlp_lower_bound_inf_);
173           z_L_vals_new[izL_new] = z_L_vals_prev[izL_prev];
174           izL_new++;
175           izL_prev++;
176         }
177         if (x_u_new[i]<nlp_upper_bound_inf_) {
178           DBG_ASSERT(x_u_prev_[i]<nlp_upper_bound_inf_);
179           z_U_vals_new[izU_new] = z_U_vals_prev[izU_prev];
180           izU_new++;
181           izU_prev++;
182         }
183       }
184       else if (x_l_prev_[i]<x_u_prev_[i]) {
185         ix_prev++;
186         izL_prev++;
187         izU_prev++;
188       }
189     }
190     DBG_ASSERT(ix_prev==prev_x->Dim());
191     DBG_ASSERT(izL_prev==prev_z_L->Dim());
192     DBG_ASSERT(izU_prev==prev_z_U->Dim());
193     DBG_ASSERT(ix_new==warm_start_iterate.x()->Dim());
194     DBG_ASSERT(izL_new==warm_start_iterate.z_L()->Dim());
195     DBG_ASSERT(izU_new==warm_start_iterate.z_U()->Dim());
196 
197     // Now copy all the values for the iterates that don't change in dimension
198     DBG_ASSERT(prev_s->Dim()==warm_start_iterate.s()->Dim());
199     DBG_ASSERT(prev_y_d->Dim()==warm_start_iterate.s()->Dim());
200     DBG_ASSERT(prev_y_d->Dim()==warm_start_iterate.y_d()->Dim());
201     for (Index i=0; i<prev_s->Dim(); i++) {
202       s_vals_new[i] = s_vals_prev[i];
203       y_d_vals_new[i] = y_d_vals_prev[i];
204     }
205     DBG_ASSERT(prev_y_c->Dim()==warm_start_iterate.y_c()->Dim());
206     for (Index i=0; i<prev_y_c->Dim(); i++) {
207       y_c_vals_new[i] = y_c_vals_prev[i];
208     }
209     DBG_ASSERT(prev_v_L->Dim()==warm_start_iterate.v_L()->Dim());
210     for (Index i=0; i<prev_v_L->Dim(); i++) {
211       v_L_vals_new[i] = v_L_vals_prev[i];
212     }
213     DBG_ASSERT(prev_v_U->Dim()==warm_start_iterate.v_U()->Dim());
214     for (Index i=0; i<prev_v_U->Dim(); i++) {
215       v_U_vals_new[i] = v_U_vals_prev[i];
216     }
217 
218     return true;
219   }
220 }
221 
222