Nonlinear Programming in IronPython QuickStart Sample

Illustrates solving nonlinear programs (optimization problems with linear or nonlinear constraints) using the NonlinearProgram and related classes in IronPython.

This sample is also available in: C#, Visual Basic, F#.

Overview

This QuickStart sample demonstrates how to solve nonlinear optimization problems with constraints using the NonlinearProgram class in Numerics.NET.

The sample presents three different approaches to solving nonlinear programming problems:

  1. A chemical equilibrium problem with linear constraints, demonstrating how to:
    • Define an objective function involving logarithms
    • Set up linear constraints using sparse matrices
    • Handle variable bounds
    • Solve a real-world optimization problem
  2. A mathematical optimization problem with nonlinear constraints, showing:
    • How to build a nonlinear program from scratch
    • Adding nonlinear constraints with manually specified gradients
    • Adding constraints where gradients are approximated numerically
    • Setting initial values for the optimization
  3. The same mathematical problem solved using automatic differentiation, illustrating:
    • How to simplify implementation by using symbolic functions
    • Letting the system automatically compute required gradients
    • Comparing results with the manual approach

Each approach is fully documented with code and explanatory comments, making it easy to understand the different ways of setting up and solving nonlinear programming problems.

The code

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 Extreme Numerics.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