Nonlinear Programming QuickStart Sample (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