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