This file illustrates the effect of nonlinear constraints, nonlinear objective functions, and specifically non-convex feasible regions in variations of the WYNDOR Glass Company example given in Hillier and Lieberman, Section 13.2.

We use the Ipopt solver backend which can handle nonlinear constraints and nonlinear objective functions.

from pyomo.environ import *
from pyomo.opt import *
#opt = solvers.SolverFactory("glpk") #Works only for LP problems
opt = solvers.SolverFactory("ipopt")
model = ConcreteModel()

About Example (4): For non-convex problems such as this, there may be several local optima and the solver may return a local optimum which is not a global optimum. The solution which gets returned depends on the initial value.

For this problem, the default initialization happens to converge to the global maximum at $x=0$ and $y=7$. However, other initial values result in convergence to the local maximum $x=4$ and $y=3$.

In the code below, we pass an explicit initial value for the search to each of the decision variables.

#### Example (0): The LP Problem ####
model.x = Var(within=NonNegativeReals)
model.y = Var(within=NonNegativeReals)

model.c1 = Constraint(expr = model.x <= 4)
model.c2 = Constraint(expr = 2*model.y <= 12)
model.c3 = Constraint(expr = 3*model.x + 2*model.y <= 18)

model.z = Objective(expr = 3*model.x + 5*model.y, sense=maximize)


#### Example (1): Nonlinear Constraint ####
# model.x = Var(within=NonNegativeReals)
# model.y = Var(within=NonNegativeReals)

# model.c1 = Constraint(expr = model.x <= 4)
# #model.c2 = Constraint(expr = 2*model.y <= 12)
# model.c3 = Constraint(expr = 9*model.x**2 + 5*model.y**2 <= 216)

# model.z = Objective(expr = 3*model.x + 5*model.y, sense=maximize)


#### Examples (2), (3): Nonlinear Objective Function ####
# model.x = Var(within=NonNegativeReals)
# model.y = Var(within=NonNegativeReals)

# model.c1 = Constraint(expr = model.x <= 4)
# model.c2 = Constraint(expr = 2*model.y <= 12)
# model.c3 = Constraint(expr = 3*model.x + 2*model.y <= 18)

# #model.z = Objective(expr = 126*model.x - 9*model.x**2 + 182*model.y - 13*model.y**2, sense=maximize)
# model.z = Objective(expr = 54*model.x - 9*model.x**2 + 78*model.y - 13*model.y**2, sense=maximize)


#### Example (4): Local vs. Global Maxima ####
# model.x = Var(within=NonNegativeReals, initialize=1)#3
# model.y = Var(within=NonNegativeReals, initialize=2)#2

# model.c1 = Constraint(expr = model.x <= 4)
# model.c2 = Constraint(expr = 2*model.y <= 14)
# model.c3 = Constraint(expr = 8*model.x - model.x**2 + 14*model.y - model.y**2 <= 49)

# model.z = Objective(expr = 3*model.x + 5*model.y, sense=maximize)

results = opt.solve(model)

model.x.get_values()

{None: 2.0000000197216092}

model.y.get_values()

{None: 6.000000059164635}

model.z.expr()

36.000000354988