1#  ___________________________________________________________________________
2#
3#  Pyomo: Python Optimization Modeling Objects
4#  Copyright 2017 National Technology and Engineering Solutions of Sandia, LLC
5#  Under the terms of Contract DE-NA0003525 with National Technology and
6#  Engineering Solutions of Sandia, LLC, the U.S. Government retains certain
7#  rights in this software.
8#  This software is distributed under the 3-clause BSD License.
9#  ___________________________________________________________________________
10
11from abc import ABCMeta, abstractmethod
12from pyomo.contrib.pynumero.interfaces import pyomo_nlp, ampl_nlp
13from pyomo.contrib.pynumero.sparse import BlockMatrix, BlockVector
14import numpy as np
15import scipy.sparse
16from pyomo.common.timing import HierarchicalTimer
17
18
19class BaseInteriorPointInterface(object, metaclass=ABCMeta):
20    @abstractmethod
21    def n_primals(self):
22        pass
23
24    @abstractmethod
25    def nnz_hessian_lag(self):
26        pass
27
28    @abstractmethod
29    def primals_lb(self):
30        pass
31
32    @abstractmethod
33    def primals_ub(self):
34        pass
35
36    @abstractmethod
37    def init_primals(self):
38        pass
39
40    @abstractmethod
41    def set_primals(self, primals):
42        pass
43
44    @abstractmethod
45    def get_primals(self):
46        pass
47
48    @abstractmethod
49    def get_obj_factor(self):
50        pass
51
52    @abstractmethod
53    def set_obj_factor(self, obj_factor):
54        pass
55
56    @abstractmethod
57    def evaluate_objective(self):
58        pass
59
60    @abstractmethod
61    def evaluate_grad_objective(self):
62        pass
63
64    @abstractmethod
65    def n_eq_constraints(self):
66        pass
67
68    @abstractmethod
69    def n_ineq_constraints(self):
70        pass
71
72    @abstractmethod
73    def nnz_jacobian_eq(self):
74        pass
75
76    @abstractmethod
77    def nnz_jacobian_ineq(self):
78        pass
79
80    @abstractmethod
81    def ineq_lb(self):
82        pass
83
84    @abstractmethod
85    def ineq_ub(self):
86        pass
87
88    @abstractmethod
89    def init_duals_eq(self):
90        pass
91
92    @abstractmethod
93    def init_duals_ineq(self):
94        pass
95
96    @abstractmethod
97    def set_duals_eq(self, duals_eq):
98        pass
99
100    @abstractmethod
101    def set_duals_ineq(self, duals_ineq):
102        pass
103
104    @abstractmethod
105    def get_duals_eq(self):
106        pass
107
108    @abstractmethod
109    def get_duals_ineq(self):
110        pass
111
112    @abstractmethod
113    def evaluate_eq_constraints(self):
114        pass
115
116    @abstractmethod
117    def evaluate_ineq_constraints(self):
118        pass
119
120    @abstractmethod
121    def evaluate_jacobian_eq(self):
122        pass
123
124    @abstractmethod
125    def evaluate_jacobian_ineq(self):
126        pass
127
128    @abstractmethod
129    def init_slacks(self):
130        pass
131
132    @abstractmethod
133    def init_duals_primals_lb(self):
134        pass
135
136    @abstractmethod
137    def init_duals_primals_ub(self):
138        pass
139
140    @abstractmethod
141    def init_duals_slacks_lb(self):
142        pass
143
144    @abstractmethod
145    def init_duals_slacks_ub(self):
146        pass
147
148    @abstractmethod
149    def set_slacks(self, slacks):
150        pass
151
152    @abstractmethod
153    def set_duals_primals_lb(self, duals):
154        pass
155
156    @abstractmethod
157    def set_duals_primals_ub(self, duals):
158        pass
159
160    @abstractmethod
161    def set_duals_slacks_lb(self, duals):
162        pass
163
164    @abstractmethod
165    def set_duals_slacks_ub(self, duals):
166        pass
167
168    @abstractmethod
169    def get_slacks(self):
170        pass
171
172    @abstractmethod
173    def get_duals_primals_lb(self):
174        pass
175
176    @abstractmethod
177    def get_duals_primals_ub(self):
178        pass
179
180    @abstractmethod
181    def get_duals_slacks_lb(self):
182        pass
183
184    @abstractmethod
185    def get_duals_slacks_ub(self):
186        pass
187
188    @abstractmethod
189    def set_barrier_parameter(self, barrier):
190        pass
191
192    @abstractmethod
193    def evaluate_primal_dual_kkt_matrix(self, timer=None):
194        pass
195
196    @abstractmethod
197    def evaluate_primal_dual_kkt_rhs(self, timer=None):
198        pass
199
200    @abstractmethod
201    def set_primal_dual_kkt_solution(self, sol):
202        pass
203
204    @abstractmethod
205    def get_delta_primals(self):
206        pass
207
208    @abstractmethod
209    def get_delta_slacks(self):
210        pass
211
212    @abstractmethod
213    def get_delta_duals_eq(self):
214        pass
215
216    @abstractmethod
217    def get_delta_duals_ineq(self):
218        pass
219
220    @abstractmethod
221    def get_delta_duals_primals_lb(self):
222        pass
223
224    @abstractmethod
225    def get_delta_duals_primals_ub(self):
226        pass
227
228    @abstractmethod
229    def get_delta_duals_slacks_lb(self):
230        pass
231
232    @abstractmethod
233    def get_delta_duals_slacks_ub(self):
234        pass
235
236    def regularize_equality_gradient(self, kkt, coef, copy_kkt=True):
237        raise RuntimeError(
238            'Equality gradient regularization is necessary but no '
239            'function has been implemented for doing so.')
240
241    def regularize_hessian(self, kkt, coef, copy_kkt=True):
242        raise RuntimeError(
243            'Hessian of Lagrangian regularization is necessary but no '
244            'function has been implemented for doing so.')
245
246
247class InteriorPointInterface(BaseInteriorPointInterface):
248    def __init__(self, pyomo_model):
249        if type(pyomo_model) is str:
250            # Assume argument is the name of an nl file
251            self._nlp = ampl_nlp.AmplNLP(pyomo_model)
252        else:
253            self._nlp = pyomo_nlp.PyomoNLP(pyomo_model)
254        self._slacks = self.init_slacks()
255
256        # set the init_duals_primals_lb/ub from ipopt_zL_out, ipopt_zU_out if available
257        # need to compress them as well and initialize the duals_primals_lb/ub
258        self._init_duals_primals_lb, self._init_duals_primals_ub =\
259            self._get_full_duals_primals_bounds()
260        self._init_duals_primals_lb[np.isneginf(self._nlp.primals_lb())] = 0
261        self._init_duals_primals_ub[np.isinf(self._nlp.primals_ub())] = 0
262        self._duals_primals_lb = self._init_duals_primals_lb.copy()
263        self._duals_primals_ub = self._init_duals_primals_ub.copy()
264
265        # set the init_duals_slacks_lb/ub from the init_duals_ineq
266        # need to be compressed and set according to their sign
267        # (-) value indicates it the upper is active, while (+) indicates
268        # that lower is active
269        self._init_duals_slacks_lb = self._nlp.init_duals_ineq().copy()
270        self._init_duals_slacks_lb[self._init_duals_slacks_lb < 0] = 0
271        self._init_duals_slacks_ub = self._nlp.init_duals_ineq().copy()
272        self._init_duals_slacks_ub[self._init_duals_slacks_ub > 0] = 0
273        self._init_duals_slacks_ub *= -1.0
274
275        self._duals_slacks_lb = self._init_duals_slacks_lb.copy()
276        self._duals_slacks_ub = self._init_duals_slacks_ub.copy()
277
278        self._delta_primals = None
279        self._delta_slacks = None
280        self._delta_duals_eq = None
281        self._delta_duals_ineq = None
282        self._barrier = None
283
284    def n_primals(self):
285        return self._nlp.n_primals()
286
287    def nnz_hessian_lag(self):
288        return self._nlp.nnz_hessian_lag()
289
290    def set_obj_factor(self, obj_factor):
291        self._nlp.set_obj_factor(obj_factor)
292
293    def get_obj_factor(self):
294        return self._nlp.get_obj_factor()
295
296    def n_eq_constraints(self):
297        return self._nlp.n_eq_constraints()
298
299    def n_ineq_constraints(self):
300        return self._nlp.n_ineq_constraints()
301
302    def nnz_jacobian_eq(self):
303        return self._nlp.nnz_jacobian_eq()
304
305    def nnz_jacobian_ineq(self):
306        return self._nlp.nnz_jacobian_ineq()
307
308    def init_primals(self):
309        primals = self._nlp.init_primals()
310        return primals
311
312    def init_slacks(self):
313        slacks = self._nlp.evaluate_ineq_constraints()
314        return slacks
315
316    def init_duals_eq(self):
317        return self._nlp.init_duals_eq()
318
319    def init_duals_ineq(self):
320        return self._nlp.init_duals_ineq()
321
322    def init_duals_primals_lb(self):
323        return self._init_duals_primals_lb
324
325    def init_duals_primals_ub(self):
326        return self._init_duals_primals_ub
327
328    def init_duals_slacks_lb(self):
329        return self._init_duals_slacks_lb
330
331    def init_duals_slacks_ub(self):
332        return self._init_duals_slacks_ub
333
334    def set_primals(self, primals):
335        self._nlp.set_primals(primals)
336
337    def set_slacks(self, slacks):
338        self._slacks = slacks
339
340    def set_duals_eq(self, duals):
341        self._nlp.set_duals_eq(duals)
342
343    def set_duals_ineq(self, duals):
344        self._nlp.set_duals_ineq(duals)
345
346    def set_duals_primals_lb(self, duals):
347        self._duals_primals_lb = duals
348
349    def set_duals_primals_ub(self, duals):
350        self._duals_primals_ub = duals
351
352    def set_duals_slacks_lb(self, duals):
353        self._duals_slacks_lb = duals
354
355    def set_duals_slacks_ub(self, duals):
356        self._duals_slacks_ub = duals
357
358    def get_primals(self):
359        return self._nlp.get_primals()
360
361    def get_slacks(self):
362        return self._slacks
363
364    def get_duals_eq(self):
365        return self._nlp.get_duals_eq()
366
367    def get_duals_ineq(self):
368        return self._nlp.get_duals_ineq()
369
370    def get_duals_primals_lb(self):
371        return self._duals_primals_lb
372
373    def get_duals_primals_ub(self):
374        return self._duals_primals_ub
375
376    def get_duals_slacks_lb(self):
377        return self._duals_slacks_lb
378
379    def get_duals_slacks_ub(self):
380        return self._duals_slacks_ub
381
382    def primals_lb(self):
383        return self._nlp.primals_lb()
384
385    def primals_ub(self):
386        return self._nlp.primals_ub()
387
388    def ineq_lb(self):
389        return self._nlp.ineq_lb()
390
391    def ineq_ub(self):
392        return self._nlp.ineq_ub()
393
394    def set_barrier_parameter(self, barrier):
395        self._barrier = barrier
396
397    def pyomo_nlp(self):
398        return self._nlp
399
400    def evaluate_primal_dual_kkt_matrix(self, timer=None):
401        if timer is None:
402            timer = HierarchicalTimer()
403        timer.start('eval hess')
404        hessian = self._nlp.evaluate_hessian_lag()
405        timer.stop('eval hess')
406        timer.start('eval jac')
407        jac_eq = self._nlp.evaluate_jacobian_eq()
408        jac_ineq = self._nlp.evaluate_jacobian_ineq()
409        timer.stop('eval jac')
410
411        duals_primals_lb = self._duals_primals_lb
412        duals_primals_ub = self._duals_primals_ub
413        duals_slacks_lb = self._duals_slacks_lb
414        duals_slacks_ub = self._duals_slacks_ub
415        primals = self._nlp.get_primals()
416
417        timer.start('hess block')
418        data = (duals_primals_lb/(primals - self._nlp.primals_lb()) +
419                duals_primals_ub/(self._nlp.primals_ub() - primals))
420        n = self._nlp.n_primals()
421        indices = np.arange(n)
422        hess_block = scipy.sparse.coo_matrix((data, (indices, indices)), shape=(n, n))
423        hess_block += hessian
424        timer.stop('hess block')
425
426        timer.start('slack block')
427        data = (duals_slacks_lb/(self._slacks - self._nlp.ineq_lb()) +
428                duals_slacks_ub/(self._nlp.ineq_ub() - self._slacks))
429        n = self._nlp.n_ineq_constraints()
430        indices = np.arange(n)
431        slack_block = scipy.sparse.coo_matrix((data, (indices, indices)), shape=(n, n))
432        timer.stop('slack block')
433
434        timer.start('set block')
435        kkt = BlockMatrix(4, 4)
436        kkt.set_block(0, 0, hess_block)
437        kkt.set_block(1, 1, slack_block)
438        kkt.set_block(2, 0, jac_eq)
439        kkt.set_block(0, 2, jac_eq.transpose())
440        kkt.set_block(3, 0, jac_ineq)
441        kkt.set_block(0, 3, jac_ineq.transpose())
442        kkt.set_block(3, 1, -scipy.sparse.identity(
443                                            self._nlp.n_ineq_constraints(),
444                                            format='coo'))
445        kkt.set_block(1, 3, -scipy.sparse.identity(
446                                            self._nlp.n_ineq_constraints(),
447                                            format='coo'))
448        timer.stop('set block')
449        return kkt
450
451    def evaluate_primal_dual_kkt_rhs(self, timer=None):
452        if timer is None:
453            timer = HierarchicalTimer()
454        timer.start('eval grad obj')
455        grad_obj = self.get_obj_factor() * self.evaluate_grad_objective()
456        timer.stop('eval grad obj')
457        timer.start('eval jac')
458        jac_eq = self._nlp.evaluate_jacobian_eq()
459        jac_ineq = self._nlp.evaluate_jacobian_ineq()
460        timer.stop('eval jac')
461        timer.start('eval cons')
462        eq_resid = self._nlp.evaluate_eq_constraints()
463        ineq_resid = self._nlp.evaluate_ineq_constraints() - self._slacks
464        timer.stop('eval cons')
465
466        timer.start('grad_lag_primals')
467        grad_lag_primals = (grad_obj +
468                            jac_eq.transpose() * self._nlp.get_duals_eq() +
469                            jac_ineq.transpose() * self._nlp.get_duals_ineq() -
470                            self._barrier / (self._nlp.get_primals() - self._nlp.primals_lb()) +
471                            self._barrier / (self._nlp.primals_ub() - self._nlp.get_primals()))
472        timer.stop('grad_lag_primals')
473
474        timer.start('grad_lag_slacks')
475        grad_lag_slacks = (-self._nlp.get_duals_ineq() -
476                           self._barrier / (self._slacks - self._nlp.ineq_lb()) +
477                           self._barrier / (self._nlp.ineq_ub() - self._slacks))
478        timer.stop('grad_lag_slacks')
479
480        rhs = BlockVector(4)
481        rhs.set_block(0, grad_lag_primals)
482        rhs.set_block(1, grad_lag_slacks)
483        rhs.set_block(2, eq_resid)
484        rhs.set_block(3, ineq_resid)
485        rhs = -rhs
486        return rhs
487
488    def set_primal_dual_kkt_solution(self, sol):
489        self._delta_primals = sol.get_block(0)
490        self._delta_slacks = sol.get_block(1)
491        self._delta_duals_eq = sol.get_block(2)
492        self._delta_duals_ineq = sol.get_block(3)
493
494    def get_delta_primals(self):
495        return self._delta_primals
496
497    def get_delta_slacks(self):
498        return self._delta_slacks
499
500    def get_delta_duals_eq(self):
501        return self._delta_duals_eq
502
503    def get_delta_duals_ineq(self):
504        return self._delta_duals_ineq
505
506    def get_delta_duals_primals_lb(self):
507        res = (((self._barrier - self._duals_primals_lb * self._delta_primals) /
508                (self._nlp.get_primals() - self._nlp.primals_lb())) -
509               self._duals_primals_lb)
510        return res
511
512    def get_delta_duals_primals_ub(self):
513        res = (((self._barrier + self._duals_primals_ub * self._delta_primals) /
514                (self._nlp.primals_ub() - self._nlp.get_primals())) -
515               self._duals_primals_ub)
516        return res
517
518    def get_delta_duals_slacks_lb(self):
519        res = (((self._barrier - self._duals_slacks_lb * self._delta_slacks) /
520                (self._slacks - self._nlp.ineq_lb())) -
521               self._duals_slacks_lb)
522        return res
523
524    def get_delta_duals_slacks_ub(self):
525        res = (((self._barrier + self._duals_slacks_ub * self._delta_slacks) /
526                (self._nlp.ineq_ub() - self._slacks)) -
527               self._duals_slacks_ub)
528        return res
529
530    def evaluate_objective(self):
531        return self._nlp.evaluate_objective()
532
533    def evaluate_eq_constraints(self):
534        return self._nlp.evaluate_eq_constraints()
535
536    def evaluate_ineq_constraints(self):
537        return self._nlp.evaluate_ineq_constraints()
538
539    def evaluate_grad_objective(self):
540        return self._nlp.evaluate_grad_objective()
541
542    def evaluate_jacobian_eq(self):
543        return self._nlp.evaluate_jacobian_eq()
544
545    def evaluate_jacobian_ineq(self):
546        return self._nlp.evaluate_jacobian_ineq()
547
548    def regularize_equality_gradient(self, kkt, coef, copy_kkt=True):
549        # Not technically regularizing the equality gradient ...
550        # Replace this with a regularize_diagonal_block function?
551        # Then call with kkt matrix and the value of the perturbation?
552
553        # Use a constant perturbation to regularize the equality constraint
554        # gradient
555        if copy_kkt:
556            kkt = kkt.copy()
557        reg_coef = coef
558        ptb = (reg_coef *
559               scipy.sparse.identity(self._nlp.n_eq_constraints(),
560                                     format='coo'))
561
562        kkt.set_block(2, 2, ptb)
563        return kkt
564
565    def regularize_hessian(self, kkt, coef, copy_kkt=True):
566        if copy_kkt:
567            kkt = kkt.copy()
568
569        hess = kkt.get_block(0, 0)
570        ptb = coef * scipy.sparse.identity(self._nlp.n_primals(), format='coo')
571        hess += ptb
572        kkt.set_block(0, 0, hess)
573        return kkt
574
575    def _get_full_duals_primals_bounds(self):
576        full_duals_primals_lb = None
577        full_duals_primals_ub = None
578        # Check in case _nlp was constructed as an AmplNLP (from an nl file)
579        if (hasattr(self._nlp, 'pyomo_model') and
580            hasattr(self._nlp, 'get_pyomo_variables')):
581            pyomo_model = self._nlp.pyomo_model()
582            pyomo_variables = self._nlp.get_pyomo_variables()
583            if hasattr(pyomo_model,'ipopt_zL_out'):
584                zL_suffix = pyomo_model.ipopt_zL_out
585                full_duals_primals_lb = np.empty(self._nlp.n_primals())
586                for i,v in enumerate(pyomo_variables):
587                    if v in zL_suffix:
588                        full_duals_primals_lb[i] = zL_suffix[v]
589
590            if hasattr(pyomo_model,'ipopt_zU_out'):
591                zU_suffix = pyomo_model.ipopt_zU_out
592                full_duals_primals_ub = np.empty(self._nlp.n_primals())
593                for i,v in enumerate(pyomo_variables):
594                    if v in zU_suffix:
595                        full_duals_primals_ub[i] = zU_suffix[v]
596
597        if full_duals_primals_lb is None:
598            full_duals_primals_lb = np.ones(self._nlp.n_primals())
599
600        if full_duals_primals_ub is None:
601            full_duals_primals_ub = np.ones(self._nlp.n_primals())
602
603        return full_duals_primals_lb, full_duals_primals_ub
604
605    def load_primals_into_pyomo_model(self):
606        if not isinstance(self._nlp, pyomo_nlp.PyomoNLP):
607            raise RuntimeError('Can only load primals into a pyomo model if a pyomo model was used in the constructor.')
608
609        pyomo_variables = self._nlp.get_pyomo_variables()
610        primals = self._nlp.get_primals()
611        for i, v in enumerate(pyomo_variables):
612            v.value = primals[i]
613
614    def pyomo_model(self):
615        return self._nlp.pyomo_model()
616
617    def get_pyomo_variables(self):
618        return self._nlp.get_pyomo_variables()
619
620    def get_pyomo_constraints(self):
621        return self._nlp.get_pyomo_constraints()
622
623    def variable_names(self):
624        return self._nlp.variable_names()
625
626    def constraint_names(self):
627        return self._nlp.constraint_names()
628
629    def get_primal_indices(self, pyomo_variables):
630        return self._nlp.get_primal_indices(pyomo_variables)
631
632    def get_constraint_indices(self, pyomo_constraints):
633        return self._nlp.get_constraint_indices(pyomo_constraints)
634