Illustrates solving nonlinear programs (optimization problems with linear or nonlinear constraints) using the NonlinearProgram and related classes in IronPython.
import numerics
from System import Array
from math import log, exp
# Vectors and matrices are in the Extreme.Mathematics namespace
from Extreme.Mathematics import *
from Extreme.Mathematics.LinearAlgebra import *
# The nonlinear programming classes reside in a namespace with
# other optimization-related classes.
from Extreme.Mathematics.Optimization import *
# Illustrates solving nonlinear programming problems
# using the classes in the Extreme.Mathematics.Optimization
# namespace of the Extreme Optimization Numerical Libraries for .NET.
# This QuickStart Sample illustrates the two ways to create a Nonlinear Program.
# The first is in terms of matrices. The coefficients
# are supplied as a matrix. The cost vector, right-hand side
# and constraints on the variables are supplied as a vector.
print "Problem with only linear constraints:"
# The variables are the concentrations of each chemical compound:
# H, H2, H2O, N, N2, NH, NO, O, O2, OH
# The objective function is the free energy, which we want to minimize:
c = Vector.Create[float](\
-6.089, -17.164, -34.054, -5.914, -24.721, \
-14.986, -24.100, -10.708, -26.662, -22.179)
objectiveFunction = lambda x: x.DotProduct(c + Vector.Log(x).AddInPlace(-log(x.Sum())))
def g1(x,y):
s = x.Sum()
(c + Vector.Log(x).AddInPlace(1.0 - log(s)) - x / s).CopyTo(y)
return y
objectiveGradient = g1
# The constraints are the mass balance relationships for each element.
# The rows correspond to the elements H, N, and O.
# The columns are the index of the variable.
# The value is the number of times the element occurs in the
# compound corresponding to the variable:
# H, H2, H2O, N, N2, NH, NO, O, O2, OH
# All this is stored in a sparse matrix, so 0 values are omitted:
A = Matrix.CreateSparse(3, 10, \
Array[int]([ 0, 0, 0, 0, 0, 1, 1, 1, 1, 2, 2, 2, 2, 2 ]), \
Array[int]([ 0, 1, 2, 5, 9, 3, 4, 5, 6, 2, 6, 7, 8, 9 ]), \
Array[float]([ 1, 2, 2, 1, 1, 1, 2, 1, 1, 1, 1, 1, 2, 1 ]))
# The right-hand sides are the atomic weights of the elements
# in the mixture.
b = Vector.Create[float](2, 1, 1)
# The number of moles for each compound must be positive.
l = Vector.CreateConstant(10, 1e-6)
u = Vector.CreateConstant(10, float.PositiveInfinity)
# We create the nonlinear program with the specified constraints:
# Because we have variable bounds, we use the constructor
# that lets us do this.
nlp1 = NonlinearProgram(objectiveFunction, objectiveGradient, A, b, b, l, u)
# We could add more (linear or nonlinear) constraints here, # but this is all we have in our problem.
nlp1.InitialGuess = Vector.CreateConstant(10, 0.1)
solution = nlp1.Solve()
print "Solution: {0:.3f}".format(solution)
print "Optimal value: {0:.5f}".format(nlp1.OptimalValue)
print "# iterations:", nlp1.SolutionReport.IterationsNeeded
print
# The second method is building the nonlinear program from scratch.
print "Problem with nonlinear constraints:"
# We start by creating a nonlinear program object. We supply
# the number of variables in the constructor.
nlp2 = NonlinearProgram(2)
nlp2.ObjectiveFunction = lambda x: exp(x[0]) * (4 * x[0] * x[0] + 2 * x[1] * x[1] + 4 * x[0] * x[1] + 2 * x[1] + 1)
def g2(x, y):
exp1 = exp(x[0])
y[0] = exp1 * (4 * x[0] * x[0] + 2 * x[1] * x[1] + 4 * x[0] * x[1] + 8 * x[0] + 6 * x[1] + 1)
y[1] = exp1 * (4 * x[0] + 4 * x[1] + 2)
return y
nlp2.ObjectiveGradient = g2
# Add constraint x0*x1 - x0 -x1 <= -1.5
def c2gradient(x, y):
y[0] = x[1] - 1
y[1] = x[0] - 1
return y
nlp2.AddNonlinearConstraint(lambda x: x[0] * x[1] - x[0] - x[1] + 1.5, \
ConstraintType.LessThanOrEqual, 0.0, c2gradient)
# Add constraint x0*x1 >= -10
# If the gradient is omitted, it is approximated using divided differences.
nlp2.AddNonlinearConstraint(lambda x: x[0] * x[1], ConstraintType.GreaterThanOrEqual, -10.0)
nlp2.InitialGuess = Vector.Create[float]( -1, 1 )
solution = nlp2.Solve()
print "Solution: {0:.6f}".format(solution)
print "Optimal value: {0:.6f}".format(nlp2.OptimalValue)
print "# iterations:", nlp2.SolutionReport.IterationsNeeded