Constraint types - Logical, general, SOS, quadratic

• Stating logic clauses that resemble SAT-type formulations (BoolVars).
• Formulating some constraints on the minimum and absolute values of linear combinations of variables (GeneralConstraints).
• Using the 'pwl' construct to formulate a piecewise linear cost function (PiecewiseLinear).
• Formulation of a small quadratic programming problem (QuadraticProgramming).
• Approximation of a nonlinear function by a special ordered set of type 2 (SpecialOrderedSets) and of a quadratic function in 2 variables by special ordered sets of type 2 (SpecialOrderedSetsQuadratic).

// (c) 2023-2024 Fair Isaac Corporation

using System;
using System.Linq;
using System.Collections.Generic;
using Optimizer.Maps;
using Optimizer.Objects;
using static Optimizer.Objects.Utils;
using static Optimizer.Objects.LinExpression;
using static Optimizer.Objects.ConstantExpression;
using Optimizer;

namespace XpressExamples
{
    /// <summary>
    /// Stating logic clauses that resemble SAT-type formulations
    /// </summary>
    public class BoolVars
    {
        public static void Main(string[] args)
        {
            Console.WriteLine("Formulating the bool vars example problem");

            const int R = 5;

            using (XpressProblem prob = new XpressProblem())
            {
                // create boolean variables and their negations
                Variable[] x = prob.AddVariables(R)
                  .WithType(ColumnType.Binary)
                  .WithName("x_{0}")
                  .WithUB(1)
                  .ToArray();
                Variable[] x_neg = prob.AddVariables(R)
                  .WithType(ColumnType.Binary)
                  .WithName("x_neg_{0}")
                  .WithUB(1)
                  .ToArray();


                // add 2 helper variables for stating logical constraints
                Variable trueVar = prob.AddVariable(0, 1, ColumnType.Binary, "TRUE");
                Variable falseVar = prob.AddVariable(0, 1, ColumnType.Binary, "FALSE");

                // fix these helper variables to TRUE and FALSE, respectively.
                prob.AddConstraint(trueVar == 1);
                prob.AddConstraint(falseVar == 0);

                // add the complement relation between each binary variable and its negation.
                prob.AddConstraints(R,
                    r => (1 * x[r] + 1 * x_neg[r] == 1));

                // add a direct equation constraint between certain variables.
                prob.AddConstraint(x[2] == x_neg[3]);

                // express some clauses on the variables
                // require that  x(0) and x_neg(4) = true
                // we use trueVar as the resultant of a logical "and" constraint
                prob.AddConstraint(trueVar.AndOf(new Variable[] { x[0], x_neg[4] }));

                // ! 'x(0) or x(2)' must be true
                prob.AddConstraint(trueVar.OrOf(new Variable[] { x[0], x[2] }));

                //  for more complicated expressions, we need auxiliary variables.
                Variable andResultant1 = prob.AddVariable(0, 1, ColumnType.Binary, "andresultant1");
                Variable andResultant2 = prob.AddVariable(0, 1, ColumnType.Binary, "andresultant2");

                // both constraints below could be formulated using andResultant[12].AndOf(...)
                // here, however, we highlight the connection to general constraints
                // the first AND is between x[0] .. x[2]
                prob.AddConstraint(andResultant1.AndOf(x.Take(3).ToArray()));
                //  the second AND is between x_neg[3] and x_neg[4]
                prob.AddConstraint(andResultant2.AndOf(x_neg.Skip(3).ToArray()));

                // add those two individual AND resultants by requiring at least one to be satisfied
                prob.AddConstraint(trueVar.OrOf(new Variable[] { andResultant1, andResultant2 }));

                // finally, add a constraint that none of x_neg[0 .. 2] should be true
                prob.AddConstraint(falseVar.OrOf(x_neg.Take(3).ToArray()));

                // set constant objective function with a maximization sense
                prob.SetObjective(Constant(0), Optimizer.ObjSense.Maximize);

                // write the problem in LP format for manual inspection
                Console.WriteLine("Writing the problem to 'BoolVars.lp'");
                prob.WriteProb("BoolVars.lp", "l");

                // Solve the problem
                Console.WriteLine("Solving the problem");
                prob.Optimize();

                // check the solution status
                Console.WriteLine("Problem finished with SolStatus {0}", prob.SolStatus);
                if (prob.SolStatus != Optimizer.SolStatus.Optimal)
                {
                    throw new Exception("Problem not solved to optimality");
                }

                // print the optimal solution of the problem to the console
                Console.WriteLine("Solution has objective value (profit) of {0}", prob.ObjVal);
                Console.WriteLine("");
                Console.WriteLine("*** Solution ***");
                double[] sol = prob.GetSolution();

                for (int r = 0; r < R; r++)
                {
                    string delim = r < R - 1 ? ", " : "\n";
                    Console.Write($"x_{r} = {x[r].GetValue(sol)}{delim}");
                }

                for (int r = 0; r < R; r++)
                {
                    string delim = r < R - 1 ? ", " : "\n";
                    Console.Write($"x_neg_{r} = {x_neg[r].GetValue(sol)}{delim}");
                }

                Console.WriteLine("");
            }
        }
    }
}

<Project Sdk="Microsoft.NET.Sdk">

  <PropertyGroup>
    <OutputType>Exe</OutputType>
    <TargetFramework>net5.0</TargetFramework>

    <IsPackable>false</IsPackable>
    <XpressExampleFiles Condition="'$(XpressExampleFiles)'==''">../../../data</XpressExampleFiles>
  </PropertyGroup>

  <ItemGroup>
    <PackageReference Include="FICO.Xpress.XPRSdn" Version="41.1.1" /> <!-- Version 41.1.1 or later -->
  </ItemGroup>

</Project>

// (c) 2023-2024 Fair Isaac Corporation

using System;
using System.Linq;
using System.Collections.Generic;
using Optimizer.Maps;
using Optimizer.Objects;
using static Optimizer.Objects.Utils;
using static Optimizer.Objects.LinExpression;
using static Optimizer.XPRSprob;
using static Optimizer.Objects.ConstantExpression;
using Optimizer;

namespace XpressExamples
{
    /// <summary>
    /// A simple example that formulates some constraints on the minimum and absolute values
    /// of linear combinations of variables.
    /// </summary>
    public class GeneralConstraints
    {
        public static void Main(string[] args)
        {
            Console.WriteLine("Formulating the general constraint example problem");

            const int R = 3;

            using (XpressProblem prob = new XpressProblem())
            {
                Variable[] x = prob.AddVariables(R)
                  .WithName("x_{0}")
                  .WithUB(20)
                  .ToArray();

                Variable y = prob.AddVariable("y");
                Variable z = prob.AddVariable("z");

                Expression objective = Sum(x);

                // We want to formulate abs(x(0)-2*x(1)) <= 10. We need to introduce
                // two auxiliary variables for the argument of the abs function
                //  and then break down the abs expression in multiple steps
                Variable diff1 = prob.AddVariable("diff1");
                Variable absOfDiff1 = prob.AddVariable("absOfDiff1");

                prob.AddConstraint(diff1 == x[0] - 2 * x[1]);
                prob.AddConstraint(absOfDiff1.AbsOf(diff1));
                prob.AddConstraint(absOfDiff1 <= 10);

                // We link a new variable to the minimum of the x(i) and
                // require this variable to be >= 5
                // Clearly, this bound constraint could also be achieved by simply setting
                // the bounds of each x variable.
                Variable minOfX = prob.AddVariable("minOfX");
                prob.AddConstraint(minOfX.MinOf(x));
                prob.AddConstraint(minOfX >= 5);

                // We link variable y to the maximum of other variables, expressions, and constants
                // y = max(x(2), 20, x(0)-z)
                Variable diff2 = prob.AddVariable("diff2");
                prob.AddConstraint(diff2 == x[0] - 1 * z);

                // the below code is equivalent to using the MaxOf function on the resultant y
                // prob.AddConstraint(y.MaxOf(new Variable[] {x[2], diff2}, 20).SetName("max_constraint"));
                prob.AddConstraint(y.MaxOf(new Variable[] { x[2], diff2 }, new double[] { 20 }, "max_constraint"));

                // set objective function with a maximization sense
                prob.SetObjective(objective, Optimizer.ObjSense.Maximize);

                // write the problem in LP format for manual inspection
                Console.WriteLine("Writing the problem to 'GeneralConstraints.lp'");
                prob.WriteProb("GeneralConstraints.lp", "l");

                // Solve the problem
                Console.WriteLine("Solving the problem");
                prob.Optimize();

                // check the solution status
                Console.WriteLine("Problem finished with SolStatus {0}", prob.SolStatus);
                if (prob.SolStatus != Optimizer.SolStatus.Optimal)
                {
                    throw new Exception("Problem not solved to optimality");
                }

                // print the optimal solution of the problem to the console
                Console.WriteLine("Solution has objective value (profit) of {0}", prob.ObjVal);
                Console.WriteLine("");
                Console.WriteLine("*** Solution ***");
                double[] sol = prob.GetSolution();

                for (int r = 0; r < R; r++)
                {
                    string delim = r < R - 1 ? ", " : "\n";
                    Console.Write($"x_{r} = {x[r].GetValue(sol)}{delim}");
                }
                Console.WriteLine($"y = {y.GetValue(sol)}, z = {z.GetValue(sol)}");
                Console.WriteLine($"ABS ( x[0] - 2*x[1] ) = {absOfDiff1.GetValue(sol)}, minOfX = {minOfX.GetValue(sol)}");

                Console.WriteLine("");
            }
        }
    }
}

<Project Sdk="Microsoft.NET.Sdk">

  <PropertyGroup>
    <OutputType>Exe</OutputType>
    <TargetFramework>net5.0</TargetFramework>

    <IsPackable>false</IsPackable>
    <XpressExampleFiles Condition="'$(XpressExampleFiles)'==''">../../../data</XpressExampleFiles>
  </PropertyGroup>

  <ItemGroup>
    <PackageReference Include="FICO.Xpress.XPRSdn" Version="41.1.1" /> <!-- Version 41.1.1 or later -->
  </ItemGroup>

</Project>

// (c) 2023-2024 Fair Isaac Corporation

using System;
using System.Linq;
using System.Collections.Generic;
using Optimizer.Maps;
using Optimizer.Objects;
using static Optimizer.Objects.Utils;
using static Optimizer.Objects.LinExpression;
using static Optimizer.XPRSprob;
using static Optimizer.Objects.ConstantExpression;
using Optimizer;

namespace XpressExamples
{
    /// <summary>
    /// An example that demonstrates how to model a piecewise linear cost function.
    /// A piecewise linear cost function f(x) assigns a different linear cost function
    /// for different intervals of the domain of its argument x.
    /// This situation occurs in real-world problems if
    /// there are price discounts available starting from a certain quantity of sold goods.
    ///
    /// - Example discussed in mipformref whitepaper -
    /// </summary>
    public class PiecewiseLinear
    {
        public static void Main(string[] args)
        {

            const int NB = 6;

            // define breakpoints as (x,y)-coordinates. Note that
            // some breakpoints share the same x-coordinate, but have
            // different values for y.
            // The optimizer handles this case correctly.
            // In practice, such noncontinuous piecewise linear functions
            // can model price discounts for example.
            PwlBreakpoint[] breakpoints = new PwlBreakpoint[NB] {
              new PwlBreakpoint(  0,   0),
              new PwlBreakpoint( 50,  50),
              new PwlBreakpoint( 50,  75),
              new PwlBreakpoint(120, 180),
              new PwlBreakpoint(120, 240),
              new PwlBreakpoint(200, 400)
            };

            Console.WriteLine("Formulating the piecewise linear example problem");

            using (XpressProblem prob = new XpressProblem())
            {
                Variable x = prob.AddVariable("x");
                Variable fx = prob.AddVariable("fx");
                Expression objective = fx;

                prob.AddConstraint(fx.PwlOf(x, breakpoints, "pwl_with_breakpoints"));

                // ! Add a lower bound on x to get a somewhat more meaningful model
                // x >= 150
                x.SetLB(150);


                // set objective function with a minimization sense
                prob.SetObjective(objective, Optimizer.ObjSense.Minimize);

                // write the problem in LP format for manual inspection
                Console.WriteLine("Writing the problem to 'PiecewiseLinear.lp'");
                prob.WriteProb("PiecewiseLinear.lp", "l");

                // Solve the problem
                Console.WriteLine("Solving the problem");
                prob.Optimize();

                // check the solution status
                Console.WriteLine("Problem finished with SolStatus {0}", prob.SolStatus);
                if (prob.SolStatus != Optimizer.SolStatus.Optimal)
                {
                    throw new Exception("Problem not solved to optimality");
                }

                // print the optimal solution of the problem to the console
                Console.WriteLine("Solution has objective value (profit) of {0}", prob.ObjVal);
                Console.WriteLine("");
                Console.WriteLine("*** Solution ***");
                double[] sol = prob.GetSolution();

                Console.WriteLine($"x = {x.GetValue(sol)}, fx = {fx.GetValue(sol)}");

                Console.WriteLine("");
            }
        }
    }
}

<Project Sdk="Microsoft.NET.Sdk">

  <PropertyGroup>
    <OutputType>Exe</OutputType>
    <TargetFramework>net5.0</TargetFramework>

    <IsPackable>false</IsPackable>
    <XpressExampleFiles Condition="'$(XpressExampleFiles)'==''">../../../data</XpressExampleFiles>
  </PropertyGroup>

  <ItemGroup>
    <PackageReference Include="FICO.Xpress.XPRSdn" Version="41.1.1" /> <!-- Version 41.1.1 or later -->
  </ItemGroup>

</Project>

// (c) 2023-2024 Fair Isaac Corporation

using Optimizer;
using Optimizer.Objects;
using static Optimizer.Objects.Utils;
using System;

namespace XpressExamples
{

    /// <summary>Small Quadratic Programming example.</summary>
    /// <remarks>
    /// <code>
    ///    minimize x1 + x1^2 +2x1x2 +2x2^2 +x4^2
    ///    s.t.
    ///      C1:  x1 +2x2 -4x4 &gt;= 0
    ///      C2: 3x1 -2x3 - x4 &lt;= 100
    ///      C3: 10 &lt;= x1 +3x2 +3x3 -2x4 &lt;= 30
    ///      0 &lt;= x1 &lt;= 20
    ///      0 &lt;= x2,x3
    ///      x4 free
    /// </code>
    /// </remarks>
    public class QuadraticProgramming
    {
        static readonly int N = 4;

        public static void Main(string[] args)
        {
            using (XpressProblem prob = new XpressProblem())
            {
                prob.callbacks.AddMessageCallback(DefaultMessageListener.Console);

                ///// VARIABLES
                Variable[] x = new Variable[N];
                x[0] = prob.AddVariable(0, 20, ColumnType.Continuous, "x1");
                x[1] = prob.AddVariable("x2");
                x[2] = prob.AddVariable("x3");
                x[3] = prob.AddVariable(Double.NegativeInfinity, Double.PositiveInfinity, ColumnType.Continuous, "x4");

                ///// OBJECTIVE
                QuadExpression obj = QuadExpression.Create();
                obj.AddTerm(x[0]);
                obj.AddTerm(x[0], x[0]);
                obj.AddTerm(2 * x[0] * x[1]);
                obj.AddTerm(2 * x[1] * x[1]);
                obj.AddTerm(x[3].Square());
                prob.SetObjective(obj, ObjSense.Minimize);

                ///// CONSTRAINTS
                prob.AddConstraint((x[0] + 2 * x[1] - 4 * x[3] >= 0).SetName("C1"));
                prob.AddConstraint((3 * x[0] - 2 * x[2] - x[3] <= 100).SetName("C2"));
                prob.AddConstraint((x[0] + 3 * x[1] + 3 * x[2] - 2 * x[3]).In(10, 30).SetName("C3"));

                ///// SOLVING + OUTPUT
                prob.WriteProb("qp.lp");
                prob.LpOptimize();

                Console.WriteLine("Problem status: " + prob.MIPStatus);
                if (prob.LPStatus != LPStatus.Optimal)
                    throw new Exception("optimization failed with status " + prob.LPStatus);

                Console.WriteLine("Objective function value: " + prob.ObjVal);
                double[] sol = prob.GetSolution();
                for (int i = 0; i < N; i++)
                    Console.Write(x[i].GetName() + ": " + x[i].GetValue(sol) + ", ");
                Console.WriteLine();
            }
        }
    }
}

<Project Sdk="Microsoft.NET.Sdk">

  <PropertyGroup>
    <OutputType>Exe</OutputType>
    <TargetFramework>net5.0</TargetFramework>

    <IsPackable>false</IsPackable>
    <XpressExampleFiles Condition="'$(XpressExampleFiles)'==''">../../../data</XpressExampleFiles>
  </PropertyGroup>

  <ItemGroup>
    <PackageReference Include="FICO.Xpress.XPRSdn" Version="41.1.1" /> <!-- Version 41.1.1 or later -->
  </ItemGroup>
  
</Project>

// (c) 2023-2024 Fair Isaac Corporation

using System;
using System.Linq;
using System.Collections.Generic;
using Optimizer.Maps;
using Optimizer.Objects;
using static Optimizer.Objects.Utils;
using static Optimizer.XPRSprob;
using static Optimizer.Objects.SOS;
using Optimizer;

namespace XpressExamples
{
    /// <summary>
    /// Approximation of a nonlinear function by a special ordered set (SOS-2).
    /// </summary>
    /// <remarks>
    /// An SOS-2 is a constraint that allows at most 2 of its variables to
    /// have a nonzero value. In addition, these variables have to be adjacent.
    ///
    /// Example discussed in mipformref whitepaper
    /// </remarks>
    public class SpecialOrderedSets
    {
        public static void Main(string[] args)
        {

            const int NB = 4; // number of breakpoints
            double[] B_X = // X coordinates of breakpoints
                new double[NB] { 1, 2.5, 4.5, 6.5 };

            double[] B_Y = // Y coordinates of breakpoints
                new double[NB] { 1.5, 6, 3.5, 2.5 };


            Console.WriteLine("Formulating the special ordered sets example problem");
            using (XpressProblem prob = new XpressProblem())
            {
                // create one w variable for each breakpoint. We express
                Variable[] w = prob.AddVariables(NB)
                  .WithName("w_{0}")
                  .WithUB(1)
                  .ToArray();

                Variable x = prob.AddVariable("x");
                Variable y = prob.AddVariable("y");

                // Define the SOS-2 with weights from B_X. This is necessary to
                // establish the ordering between the w variables.
                prob.AddConstraint(SOS.Sos2(w, B_X, "SOS_2"));

                // We use the w variables to express a convex combination.
                // In combination with the above SOS-2 condition,
                // X and Y are represented as a convex combination of 2 adjacent
                // breakpoints.
                prob.AddConstraint(Sum(w) == 1);

                // The below constraints express the actual locations of X and Y
                // in the plane as a convex combination of the breakpoints, subject
                // to the assignment found for the w variables.
                prob.AddConstraint(x == ScalarProduct(w, B_X));
                prob.AddConstraint(y == ScalarProduct(w, B_Y));

                // set lower and upper bounds on x
                x.SetLB(2); x.SetUB(6);

                // set objective function with a maximization sense
                prob.SetObjective(y, Optimizer.ObjSense.Minimize);

                // write the problem in LP format for manual inspection
                Console.WriteLine("Writing the problem to 'SpecialOrderedSets.lp'");
                prob.WriteProb("SpecialOrderedSets.lp", "l");

                // Solve the problem
                Console.WriteLine("Solving the problem");
                prob.Optimize();

                // check the solution status
                Console.WriteLine("Problem finished with SolStatus {0}", prob.SolStatus);
                if (prob.SolStatus != Optimizer.SolStatus.Optimal)
                {
                    throw new Exception("Problem not solved to optimality");
                }

                // print the optimal solution of the problem to the console
                Console.WriteLine("Solution has objective value (profit) of {0}\n", prob.ObjVal);
                Console.WriteLine("*** Solution ***");
                double[] sol = prob.GetSolution();

                for (int b = 0; b < NB; b++)
                {
                    string delim = b < NB - 1 ? ", " : "\n";
                    Console.Write($"w_{b} = {w[b].GetValue(sol)}{delim}");
                }

                Console.WriteLine($"x = {x.GetValue(sol)}, y = {y.GetValue(sol)}");
            }
        }
    }
}

<Project Sdk="Microsoft.NET.Sdk">

  <PropertyGroup>
    <OutputType>Exe</OutputType>
    <TargetFramework>net5.0</TargetFramework>

    <IsPackable>false</IsPackable>
    <XpressExampleFiles Condition="'$(XpressExampleFiles)'==''">../../../data</XpressExampleFiles>
  </PropertyGroup>

  <ItemGroup>
    <PackageReference Include="FICO.Xpress.XPRSdn" Version="41.1.1" /> <!-- Version 41.1.1 or later -->
  </ItemGroup>

</Project>

// (c) 2023-2024 Fair Isaac Corporation

using System;
using System.Linq;
using System.Collections.Generic;
using Optimizer.Maps;
using Optimizer.Objects;
using static Optimizer.Objects.Utils;
using static Optimizer.XPRSprob;
using static Optimizer.Objects.ConstantExpression;
using static Optimizer.Objects.SOS;
using Optimizer;

namespace XpressExamples
{
    /// <summary>
    /// Approximation of a quadratic function in 2 variables by special ordered sets (SOS-2).
    /// </summary>
    /// <remarks>
    /// An SOS-2 is a constraint that allows at most 2 of its variables to
    /// have a nonzero value. In addition, these variables have to be adjacent.
    ///
    /// Example discussed in mipformref whitepaper.
    /// </remarks>
    public class SpecialOrderedSetsQuadratic
    {
        public static void Main(string[] args)
        {

            const int NX = 10; // number of breakpoints on the X-axis
            const int NY = 10; // number of breakpoints on the Y-axis
            double[] X = // X coordinates of grid points
                new double[NX];

            double[] Y = // Y coordinates of breakpoints
                new double[NY];

            double[,] F_XY = // two dimensional array of function values on the grid points
                new double[NX, NY];

            // assign the toy data
            for (int i = 0; i < NX; i++) X[i] = i + 1;
            for (int i = 0; i < NY; i++) Y[i] = i + 1;
            for (int i = 0; i < NX; i++) for (int j = 0; j < NY; j++) F_XY[i, j] = (X[i] - 5) * (Y[j] - 5);

            Console.WriteLine("Formulating the special ordered sets quadratic example problem");
            using (XpressProblem prob = new XpressProblem())
            {
                // create one w variable for each X breakpoint. We express
                Variable[] wx = prob.AddVariables(NX)
                  .WithName("wx_{0}")
                  .WithUB(1)  // this upper bound i redundant because of the convex combination constraint on the sum of the wx
                  .ToArray();
                // create one w variable for each Y breakpoint. We express
                Variable[] wy = prob.AddVariables(NY)
                  .WithName("wy_{0}")
                  .WithUB(1)  // this upper bound i redundant because of the convex combination constraint on the sum of the wy
                  .ToArray();

                // create a two-dimensional array of w variable for each grid point. We express
                Variable[,] wxy = prob.AddVariables(NX, NY)
                  .WithName("wxy_{0}_{1}")
                  .WithUB(1)  // this upper bound i redundant because of the convex combination constraint on the sum of the wy
                  .ToArray();

                Variable x = prob.AddVariable("x");
                Variable y = prob.AddVariable("y");
                Variable fxy = prob.AddVariable("fxy");

                // make fxy a free variable
                fxy.SetLB(double.NegativeInfinity);

                // Define the SOS-2 constraints with weights from X and Y. This
                // is necessary to establish the ordering between variables in wx
                // and in wy.
                prob.AddConstraint(SOS.Sos2(wx, X, "SOS_2_X"));
                prob.AddConstraint(SOS.Sos2(wy, Y, "SOS_2_Y"));

                prob.AddConstraint(Sum(wx) == 1);
                prob.AddConstraint(Sum(wy) == 1);

                // link the wxy variables to their 1-dimensional colleagues
                prob.AddConstraints(NX,
                  i => (wx[i] == Sum(NY, j => wxy[i, j]))
                );
                prob.AddConstraints(NY,
                  j => (wy[j] == Sum(NX, i => wxy[i, j]))
                );


                // now express the actual x, y, and f(x,y) coordinates
                prob.AddConstraint(x == ScalarProduct(wx, X));
                prob.AddConstraint(y == ScalarProduct(wy, Y));
                prob.AddConstraint(fxy == Sum(NX, i => Sum(NY, j => wxy[i, j] * F_XY[i, j])));

                // set lower and upper bounds on x and y
                x.SetLB(2); x.SetUB(10);
                y.SetLB(2); y.SetUB(10);

                // set objective function with a minimization sense
                prob.SetObjective(fxy, Optimizer.ObjSense.Minimize);

                // write the problem in LP format for manual inspection
                Console.WriteLine("Writing the problem to 'SpecialOrderedSetsQuadratic.lp'");
                prob.WriteProb("SpecialOrderedSetsQuadratic.lp", "l");

                // Solve the problem
                Console.WriteLine("Solving the problem");
                prob.Optimize();

                // check the solution status
                Console.WriteLine("Problem finished with SolStatus {0}", prob.SolStatus);
                if (prob.SolStatus != Optimizer.SolStatus.Optimal)
                {
                    throw new Exception("Problem not solved to optimality");
                }

                // print the optimal solution of the problem to the console
                Console.WriteLine("Solution has objective value (profit) of {0}\n", prob.ObjVal);
                Console.WriteLine("*** Solution ***");
                double[] sol = prob.GetSolution();

                for (int i = 0; i < NX; i++)
                {
                    string delim = i < NX - 1 ? ", " : "\n";
                    Console.Write($"wx_{i} = {wx[i].GetValue(sol)}{delim}");
                }
                for (int j = 0; j < NY; j++)
                {
                    string delim = j < NY - 1 ? ", " : "\n";
                    Console.Write($"wy_{j} = {wy[j].GetValue(sol)}{delim}");
                }

                Console.WriteLine($"x = {x.GetValue(sol)}, y = {y.GetValue(sol)}");
            }
        }
    }
}