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 pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP
12from pyomo.contrib.pynumero.interfaces.utils import (build_bounds_mask,
13                                                     build_compression_matrix,
14                                                     full_to_compressed)
15import pyomo.environ as pyo
16import numpy as np
17
18
19def create_basic_model():
20
21    m = pyo.ConcreteModel()
22    m.x = pyo.Var([1, 2, 3], domain=pyo.Reals)
23    for i in range(1, 4):
24        m.x[i].value = i
25    m.c1 = pyo.Constraint(expr=m.x[1] ** 2 - m.x[2] - 1 == 0)
26    m.c2 = pyo.Constraint(expr=m.x[1] - m.x[3] - 0.5 == 0)
27    m.d1 = pyo.Constraint(expr=m.x[1] + m.x[2] <= 100.0)
28    m.d2 = pyo.Constraint(expr=m.x[2] + m.x[3] >= -100.0)
29    m.d3 = pyo.Constraint(expr=m.x[2] + m.x[3] + m.x[1] >= -500.0)
30    m.x[2].setlb(0.0)
31    m.x[3].setlb(0.0)
32    m.x[2].setub(100.0)
33    m.obj = pyo.Objective(expr=m.x[2]**2)
34    return m
35
36
37def main():
38    model = create_basic_model()
39    solver = pyo.SolverFactory('ipopt')
40    solver.solve(model, tee=True)
41
42    # build nlp initialized at the solution
43    nlp = PyomoNLP(model)
44
45    # get initial point
46    print(nlp.primals_names())
47    x0 = nlp.get_primals()
48
49    # vectors of lower and upper bounds
50    xl = nlp.primals_lb()
51    xu = nlp.primals_ub()
52
53    # demonstrate use of compression from full set of bounds
54    # to only finite bounds using masks
55    xlb_mask = build_bounds_mask(xl)
56    xub_mask = build_bounds_mask(xu)
57    # get the compressed vector
58    compressed_xl = full_to_compressed(xl, xlb_mask)
59    compressed_xu = full_to_compressed(xu, xub_mask)
60    # we can also build compression matrices
61    Cx_xl = build_compression_matrix(xlb_mask)
62    Cx_xu = build_compression_matrix(xub_mask)
63
64    # lower and upper bounds residual
65    res_xl = Cx_xl * x0 - compressed_xl
66    res_xu = compressed_xu - Cx_xu * x0
67    print("Residuals lower bounds x-xl:", res_xl)
68    print("Residuals upper bounds xu-x:", res_xu)
69
70    # set the value of the primals (we can skip the duals)
71    # here we set them to the initial values, but we could
72    # set them to anything
73    nlp.set_primals(x0)
74
75    # evaluate residual of equality constraints
76    print(nlp.constraint_names())
77    res_eq = nlp.evaluate_eq_constraints()
78    print("Residuals of equality constraints:", res_eq)
79
80    # evaluate residual of inequality constraints
81    res_ineq = nlp.evaluate_ineq_constraints()
82
83    # demonstrate the use of compression from full set of
84    # lower and upper bounds on the inequality constraints
85    # to only the finite values using masks
86    ineqlb_mask = build_bounds_mask(nlp.ineq_lb())
87    inequb_mask = build_bounds_mask(nlp.ineq_ub())
88    # get the compressed vector
89    compressed_ineq_lb = full_to_compressed(nlp.ineq_lb(), ineqlb_mask)
90    compressed_ineq_ub = full_to_compressed(nlp.ineq_ub(), inequb_mask)
91    # we can also build compression matrices
92    Cineq_ineqlb = build_compression_matrix(ineqlb_mask)
93    Cineq_inequb = build_compression_matrix(inequb_mask)
94
95    # lower and upper inequalities residual
96    res_ineq_lb = Cineq_ineqlb * res_ineq - compressed_ineq_lb
97    res_ineq_ub = compressed_ineq_ub - Cineq_inequb*res_ineq
98    print("Residuals of inequality constraints lower bounds:", res_ineq_lb)
99    print("Residuals of inequality constraints upper bounds:", res_ineq_ub)
100
101    feasible = False
102    if np.all(res_xl >= 0) and np.all(res_xu >= 0) \
103        and np.all(res_ineq_lb >= 0) and np.all(res_ineq_ub >= 0) and \
104        np.allclose(res_eq, np.zeros(nlp.n_eq_constraints()), atol=1e-5):
105        feasible = True
106
107    print("Is x0 feasible:", feasible)
108
109    return feasible
110
111
112if __name__ == '__main__':
113    main()
114