Large-scale nonlinear optimization for F#

This blog shows how to use Ipopt (https://projects.coin-or.org/Ipopt) in F#.
Ipopt (Interior Point OPTimizer) is a software package for large-scale nonlinear optimization.
It is designed to find (local) solutions of mathematical optimization problems of the from
   min f(x)
x in R^n
s.t.    g_L <= g(x) <= g_U
        x_L <= x <= x_U

where f(x): R^n --> R is the objective function, and g(x): R^n --> R^m are the constraint functions. The vectors g_L and g_U denote the lower and upper bounds on the constraints, and the vectors x_L and x_U are the bounds on the variables x. The functions f(x) and g(x) can be nonlinear and nonconvex, but should be twice continuously differentiable. Note that equality constraints can be formulated in the above formulation by setting the corresponding components of g_L and g_U to the same value.

The binaries can be downloaded from http://www.coin-or.org/download/binary/Ipopt/. I am using the library Ipopt38.dll.
The interface for .NET can be downloaded from http://code.google.com/p/csipopt/.
Following the instruction on README.txt, you can create a class library Cureos.Numeric.dll by running the batch script build_library.bat.
Based on the C# example, provided by csipopt, I created the following F# code. I created a type, following the example, but it is not really necessary.

#I @"C:\Users\caxelrud\Documents\Visual Studio 2010\Projects\hs071_fsx"
#r "Cureos.Numerics.dll"
//Need Ipopt38.dll in path

open System
open System.IO
open Cureos.Numerics
open System.Runtime.InteropServices

type HS071=
    val public n : int
    val public m : int
    val public nele_jac : int
    val public nele_hess : int
    val public x_L : double[]
    val public x_U : double[]
    val public g_L : double[]
    val public g_U : double[]
    //set the number of variables and allocate space for the bounds
    //set the values for the variable bounds using n, x_L, x_U
    //Set the number of constraints and allocate space for the bounds using m, g_L, g_U
    //Set the number of nonzeros in the Jacobian of the constraints using nele_jac
    //Set the number of nonzeros in the Hessian of the Lagrangian (lower or upper triangual part only) using nele_hess

    new ()= {n=4;m=2;nele_jac = 8;nele_hess = 10;x_L =[1.0; 1.0; 1.0; 1.0 ];x_U = [5.0; 5.0; 5.0; 5.0];g_L=[25.0; 40.0];g_U = [Ipopt.PositiveInfinity; 40.0]}
 
    member public this.eval_f (n:int, x:double[], new_x:bool, [] obj_value:byref)=
        obj_value <- x.[0]*x.[3]*(x.[0]+x.[1]+x.[2])+x.[2]
        true

    member public this.eval_grad_f (n:int, x:double[], new_x:bool, grad_f:double[])=
        grad_f.[0] <- x.[0] * x.[3] + x.[3] * (x.[0] + x.[1] + x.[2])
        grad_f.[1] <- x.[0] * x.[3]
        grad_f.[2] <- x.[0] * x.[3] + 1.0
        grad_f.[3] <- x.[0] * (x.[0] + x.[1] + x.[2])
        true

    member public this.eval_g (n:int, x:double[], new_x:bool, m:int, g:double[])=

        g.[0] <- x.[0] * x.[1] * x.[2] * x.[3];
        g.[1] <- x.[0] * x.[0] + x.[1] * x.[1] + x.[2] * x.[2] + x.[3] * x.[3];
        true;

    member public this.eval_jac_g (n:int, x:double[], new_x:bool, m:int, nele_jac:int, iRow:int[], jCol:int[], values:double[])=
        if values = null then
            // set the structure of the jacobian
           // this particular jacobian is dense
            iRow.[0] <- 0
            jCol.[0] <- 0
            iRow.[1] <- 0
            jCol.[1] <- 1
            iRow.[2] <- 0
            jCol.[2] <- 2
            iRow.[3] <- 0
            jCol.[3] <- 3
            iRow.[4] <- 1
            jCol.[4] <- 0
            iRow.[5] <- 1
           jCol.[5] <- 1
           iRow.[6] <- 1
           jCol.[6] <- 2
           iRow.[7] <- 1
           jCol.[7] <- 3
       else
          // return the values of the jacobian of the constraints
          values.[0] <- x.[1] * x.[2] * x.[3] // 0,0
          values.[1] <- x.[0] * x.[2] * x.[3] // 0,1
          values.[2] <- x.[0] * x.[1] * x.[3] // 0,2
          values.[3] <- x.[0] * x.[1] * x.[2] // 0,3
          values.[4] <- 2.0 * x.[0] // 1,0
          values.[5] <- 2.0 * x.[1] // 1,1
          values.[6] <- 2.0 * x.[2] // 1,2
          values.[7] <- 2.0 * x.[3] // 1,3
       true

    member public this.eval_h (n:int, x:double[], new_x:bool, obj_factor:double,m:int, lambda:double[],  new_lambda:bool,nele_hess:int, iRow:int[], jCol:int[],values:double[]) =

        if values = null then
            // set the Hessian structure. This is a symmetric matrix, fill the lower left
            // triangle only.
            // the hessian for this problem is actually dense
            let mutable idx= 0; // nonzero element counter
            for row in 0..3 do
                 for col in 0..row do
                      iRow.[idx] <- row
                      jCol.[idx] <- col
                      idx<-idx+1
        else
            // return the values. This is a symmetric matrix, fill the lower left
            // triangle only
            // fill the objective portion
            values.[0] <- obj_factor * (2.0 * x.[3]); // 0,0
            values.[1] <- obj_factor * (x.[3]); // 1,0
            values.[2] <- 0.0; // 1,1
            values.[3] <- obj_factor * (x.[3]); // 2,0
            values.[4] <- 0.0; // 2,1
            values.[5] <- 0.0; // 2,2
            values.[6] <- obj_factor * (2.0 * x.[0] + x.[1] + x.[2]); // 3,0
            values.[7] <- obj_factor * (x.[0]); // 3,1
            values.[8] <- obj_factor * (x.[0]); // 3,2
            values.[9] <- 0.0; // 3,3
            // add the portion for the first constraint
            values.[1] <- values.[1]+ lambda.[0] * (x.[2] * x.[3]); // 1,0 */
            values.[3] <- values.[3]+ lambda.[0] * (x.[1] * x.[3]); // 2,0 */
            values.[4] <- values.[4]+ lambda.[0] * (x.[0] * x.[3]); // 2,1 */
            values.[6] <- values.[6]+ lambda.[0] * (x.[1] * x.[2]); // 3,0 */
            values.[7] <- values.[7]+ lambda.[0] * (x.[0] * x.[2]); // 3,1 */
            values.[8] <- values.[8]+ lambda.[0] * (x.[0] * x.[1]); // 3,2 */
            // add the portion for the second constraint
            values.[0] <- values.[0]+ lambda.[1] * 2.0; // 0,0 */
            values.[2] <- values.[2]+ lambda.[1] * 2.0; // 1,1 */
            values.[5] <- values.[5]+ lambda.[1] * 2.0; // 2,2 */
            values.[9] <- values.[9]+ lambda.[1] * 2.0; // 3,3 */
        true

// create the IpoptProblem

let p = new HS071()
// allocate space for the initial point and set the values
let x = [1.0; 5.0; 5.0; 1.0 ]
let mutable obj=0.0
let d_eval_g= new EvaluateConstraintsDelegate(fun (n:int) (x:double[]) (new_x:bool) (m:int) (g:double[])-> p.eval_g(n,x,new_x,m,g))
let d_eval_f= new EvaluateObjectiveDelegate(fun n x new_x obj_value-> p.eval_f(n,x,new_x,&obj_value))
let d_eval_grad_f= new EvaluateObjectiveGradientDelegate(fun n x new_x grad_f-> p.eval_grad_f(n,x,new_x,grad_f))
let d_eval_jac_g= new EvaluateJacobianDelegate(fun n x new_x m nele_jac iRow jCol values-> p.eval_jac_g(n,x,new_x,m,nele_jac,iRow,jCol,values))
let d_eval_h= new EvaluateHessianDelegate(fun n x new_x obj_factor m lambda new_Lambda nele_hess iRow jCol values-> p.eval_h(n,x,new_x,obj_factor,m,lambda,new_Lambda,nele_hess,iRow,jCol,values))

let problem = new Ipopt(p.n, p.x_L, p.x_U, p.m, p.g_L, p.g_U, p.nele_jac, p.nele_hess, d_eval_f, d_eval_g, d_eval_grad_f, d_eval_jac_g, d_eval_h)

// Set some options. Note the following ones are only examples, they might not be suitable for your problem.
problem.AddOption("tol", 1e-7) >ignore
problem.AddOption("mu_strategy", "adaptive") >ignore
problem.AddOption("output_file", "hs071_fs.txt") >ignore

// solve the problem
let status=problem.SolveProblem(x, &obj, null, null, null, null);
printfn "Optimization return status:%A" status
    for i in 0..3 do
        printfn "x[%d]=%f" i x.[i]
(*
Optimization return status:Solve_Succeeded
x[0]=1.000000
x[1]=4.743000
x[2]=3.821150
x[3]=1.379408
*)