FICO
FICO Xpress Optimization Examples Repository
FICO Optimization Community FICO Xpress Optimization Home
Back to examples browserPrevious example

Maximizing the area of a polygon using a 'multimapdelta' userfunction

Description
Demonstrates how to solve a nonlinear problem in C#

Further explanation of this example: 'Xpress NonLinear Reference Manual'

Polygon_dnet_multimapdelta.zip[download all files]

Source Files
By clicking on a file name, a preview is opened at the bottom of this page.
PolygonMultiMapDelta.cs[download]
PolygonMultiMapDelta.csproj[download]





PolygonMultiMapDelta.cs

using System;
using Optimizer;
using System.Collections.Generic;

namespace Examples {
    /// <summary>
    /// Code example that uses a user function of type "multimapdelta".
    /// </summary>
    /// <remarks>
    /// <code>
    ///     Xpress Optimizer Examples
    ///     =========================
    ///
    ///     Maximize the area of polygon of N vertices and diameter of 1
    ///     The position of vertices is indicated as (rho,theta) coordinates
    ///     where rho denotes the distance to the base point
    ///     (vertex with number N) and theta the angle from the x-axis.
    ///
    ///     (c) 2021-2024 Fair Isaac Corporation
    /// </code>
    ///
    ///     Polygon example: maximise the area of an N sided polygon
    ///
    ///     *** Demonstrating using a multi-map (R^2->R^2) userfunction calculating its own derivatives ***
    ///
    /// <code>
    ///     Variables:
    ///
    ///     rho   : 0..N-1   ! Distance of vertex from the base point
    ///     theta : 0..N-1   ! Angle from x-axis
    ///
    ///     Objective:
    ///         (sum (i in 1..N-2) (rho(i)*rho(i-1)*sin(theta(i)-theta(i-1)))) * 0.5
    ///
    ///     Constraints:
    ///         Vertices in increasing degree order:
    ///         theta(i) >= theta(i-1) +.0001  : i = 1..N-2
    ///         Boundary conditions:
    ///         theta(N-1) <= Pi
    ///         0.1 <= rho(i) <= 1   : i = 0..N-2
    ///         Third side of all triangles <= 1
    ///         rho(i)^2 + rho(j)^2 - rho(i)*rho(j)*2*cos(theta(j)-theta(i)) <= 1    : i in 0..N-3, j in i..N-2
    /// </code>
    /// </remarks>
    public class PolygonMultiMapDelta {

        /// <summary>
        /// User function that maps two doubles to a double.
        /// </summary>
        /// <remarks>
        /// Computes <c>sin(x[0] - x[1])</c>.
        /// It also fills in the partial derivatives if requested.
        /// </remarks>
        /// <param name="value">The point at which to evaluate.</param>
        /// <param name="delta">Non-zero whether partial derivatives are requested.</param>
        /// <paran name="partials">Where to return the partial derivatives of the function.</paran>
        /// <returns>the computed values of the function.</returns>
        private static double[] MyTrigonometrics(double[] value, double[] deltas, double[] partials)
        {
            int nInput = 2;
            // Assuming f:R^k->R^l, there will be a total of k*l derivatives,
            // which must be written to the partials argument as:
            //   diff(f1(x), x1), diff(f1(x), x2) ... diff(f1(x), xk)
            //   diff(f2(x), x1), diff(f2(x), x2) ... diff(f1(x), xk)
            //   ...
            //   diff(fl(x), x1), diff(fl(x), x2) ... diff(fl(x), xk)

            if (deltas != null) // Delta may be used as a suggestion for a finite difference step size
                                // however it also indicates if a partial is requested, saving on effort in case only an evaluation is needed
            {
                if (deltas[0] != 0.0)
                {
                    partials[0 + 0] = Math.Cos(value[0] - value[1]);
                    partials[nInput + 0] = -Math.Sin(value[0] - value[1]);
                }
                if (deltas[1] != 0.0)
                {
                    partials[0 + 1] = -Math.Cos(value[0] - value[1]);
                    partials[nInput + 1] = Math.Sin(value[0] - value[1]);
                }
            }

            return new double[] {
                Math.Sin(value[0] - value[1]),
                Math.Cos(value[0] - value[1])
            };
        }

        public static void Main(String[] args) {
            using (XPRSprob prob = new XPRSprob(null)) {
                // Number of sides of the Polygon
                int nSide = 5;

                // Theta
                int[] theta = prob.VarArray('C', nSide - 1, 0.0, Math.PI,
                                            i => "THETA" + (i + 1));
                // Rho
                int[] rho = prob.VarArray('C', nSide - 1, 0.01, 1.0,
                                          i => "RHO" + (i + 1));

                // Add the user function
                XPRSprob.MultiMapDeltaFunction trig = prob.NlpAddUserFunction("myTrigonometrics", 2, 2, 0, MyTrigonometrics);

                // Objective function. We build the objective function as
                // a formula in infix notation. See below for submitting a
                // formula as string.
                List<int> tok = new List<int>();
                List<double> val = new List<double>();
                for (int i = 1; i < nSide - 1; ++i) {
                    if ( tok.Count > 0 ) {
                        tok.Add(Constants.TOK_OP);
                        val.Add(Constants.OP_PLUS);
                    }
                    tok.Add(Constants.TOK_COL);             // RHO(i)
                    val.Add(rho[i]);
                    tok.Add(Constants.TOK_OP);              // *
                    val.Add(Constants.OP_MULTIPLY);
                    tok.Add(Constants.TOK_COL);             // RHO(i-1)
                    val.Add(rho[i-1]);
                    tok.Add(Constants.TOK_OP);              // *
                    val.Add(Constants.OP_MULTIPLY);
                    tok.Add(Constants.TOK_FUN);             // myTrigonometrics
                    val.Add(trig.GetId());
                    tok.Add(Constants.TOK_LB);              // (
                    val.Add(Constants.TOK_LB);
                    tok.Add(Constants.TOK_COL);             // THETA(i)
                    val.Add(theta[i]);
                    tok.Add(Constants.TOK_DEL);             // -
                    val.Add(Constants.DEL_COMMA);
                    tok.Add(Constants.TOK_COL);             // THETA(i-1)
                    val.Add(theta[i-1]);
                    tok.Add(Constants.TOK_DEL);             // :
                    val.Add(Constants.DEL_COLON);
                    tok.Add(Constants.TOK_CON);             // 1
                    val.Add(1);
                    tok.Add(Constants.TOK_RB);              // )
                    val.Add(Constants.TOK_RB);
                    tok.Add(Constants.TOK_OP);              // *
                    val.Add(Constants.OP_MULTIPLY);
                    tok.Add(Constants.TOK_CON);             // 0.5
                    val.Add(0.5);
                }
                tok.Add(Constants.TOK_EOF);
                val.Add(0.0);
                // Since nonlinear objectives cannot be directly expressed in Xpress, we maximize a free
                // variable objx and constrain this variable to be equal to the nonlinear objective.
                int objx = prob.AddCol(1.0, XPRS.MINUSINFINITY, XPRS.PLUSINFINITY);
                int objeq = prob.AddRow(new int[]{objx}, new double[]{-1.0}, 'E', 0.0);
                prob.NlpChgFormula(objeq, 0,
                                   tok.ToArray(),
                                   val.ToArray());
                prob.ChgObjSense(ObjSense.Maximize);

                // Vertices in increasing degree order
                for (int i = 0; i < nSide - 2; ++i)
                    prob.AddRow(new int[]{ theta[i], theta[i+1] },
                                new double[]{ -1.0, 1.0 },
                                'G',
                                0.001);

                // Third side of all triangles <= 1
                for (int i = 1; i < nSide - 1; i++) {
                    for (int j = i + 1; j < nSide; j++) {
                        int row = prob.AddRow(new int[0], new double[0], 'L', 1.0);
                        prob.NlpChgFormulaStr(row,
                            String.Format("RHO{0:d} ^ 2 + RHO{1:d} ^ 2 - RHO{2:d} * RHO{3:d} * 2 * myTrigonometrics ( THETA{4:d} , THETA{5:d} : 2 )",
                            i, j, i, j, j, i));
                    }
                }

                int solvestatus = -1;
                int solstatus = -1;
                prob.Optimize("s", out solvestatus, out solstatus);
                System.Console.WriteLine("solvestatus: " + solvestatus);
                System.Console.WriteLine("solstatus:   " + solstatus);
                System.Console.WriteLine("Solution value: " + prob.ObjVal);

                double[] x = new double[prob.Cols];
                int status = -1;
                prob.GetSolution(out status, x, 0, x.Length - 1);
                for (int i = 0; i < rho.Length; ++i)
                    System.Console.WriteLine("RHO" + (i+1) + " = " + x[rho[i]]);
                for (int i = 0; i < theta.Length; ++i)
                    System.Console.WriteLine("THETA" + (i + 1) + " = " + x[theta[i]]);
            }
        }
    }
}

Back to examples browserPrevious example