16 List of examples

List of examples shipped in the distribution of Fusion API for .NET:

Table 16.1 List of distributed examples

File

Description

TrafficNetworkModel.cs

Demonstrates a traffic network problem as a conic quadratic problem (CQO)

alan.cs

A portfolio choice model alan.gms from the GAMS model library

baker.cs

A small bakery revenue maximization linear problem

breaksolver.cs

Shows how to break a long-running task

callback.cs

An example of data/progress callback

ceo1.cs

A simple conic exponential problem

cqo1.cs

A simple conic quadratic problem

diet.cs

Solving Stigler’s Nutrition model diet from the GAMS model library

djc1.cs

A simple problem with disjunctive constraints (DJC)

duality.cs

Shows how to access the dual solution

elastic.cs

Linear regression with elastic net. Demonstrates model parametrization.

facility_location.cs

Demonstrates a small one-facility location problem (CQO)

gp1.cs

A simple geometric program (GP) in conic form

lo1.cs

A simple linear problem

logistic.cs

Implements logistic regression and simple log-sum-exp (CEO)

lownerjohn_ellipsoid.cs

Computes the Lowner-John inner and outer ellipsoidal approximations of a polytope (SDO, Power Cone)

lpt.cs

Demonstrates how to solve the multi-processor scheduling problem and input an integer feasible point (MIP)

mico1.cs

A simple mixed-integer conic problem

milo1.cs

A simple mixed-integer linear problem

mioinitsol.cs

A simple mixed-integer linear problem with an initial guess

nearestcorr.cs

Solves the nearest correlation matrix problem (SDO, CQO)

opt_server_sync.cs

Uses MOSEK OptServer to solve an optimization problem synchronously

parallel.cs

Demonstrates parallel optimization using a batch method in MOSEK

parameters.cs

Shows how to set optimizer parameters and read information items

portfolio_1_basic.cs

Portfolio optimization - basic Markowitz model

portfolio_2_frontier.cs

Portfolio optimization - efficient frontier

portfolio_3_impact.cs

Portfolio optimization - market impact costs

portfolio_4_transcost.cs

Portfolio optimization - transaction costs

portfolio_5_card.cs

Portfolio optimization - cardinality constraints

portfolio_6_factor.cs

Portfolio optimization - factor model

pow1.cs

A simple power cone problem

primal_svm.cs

Implements a simple soft-margin Support Vector Machine (CQO)

qcqp_sdo_relaxation.cs

Demonstrate how to use SDP to solve convex relaxation of a mixed-integer QCQO problem

reoptimization.cs

Demonstrate how to modify and re-optimize a linear problem

response.cs

Demonstrates proper response handling

sdo1.cs

A simple semidefinite problem with one matrix variable and a quadratic cone

sdo2.cs

A simple semidefinite problem with two matrix variables

sdo3.cs

A simple semidefinite problem with many matrix variables of the same dimension

sospoly.cs

Models the cone of nonnegative polynomials and nonnegative trigonometric polynomials using Nesterov’s framework

sudoku.cs

A SUDOKU solver (MIP)

total_variation.cs

Demonstrates how to solve a total variation problem (CQO)

tsp.cs

Solves a simple Travelling Salesman Problem and shows how to add constraints to a model and re-optimize (MIP)

Additional examples can be found on the MOSEK website and in other MOSEK publications.

TrafficNetworkModel.cs

Listing 16.1 TrafficNetworkModel.cs Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      TrafficNetworkModel.cs
//
// Purpose:   Demonstrates a traffix network problem as a conic quadratic problem.
//
// Source:    Robert Fourer, "Convexity Checking in Large-Scale Optimization",
//            OR 53 --- Nottingham 6-8 September 2011.
//
// The problem:
//            Given a directed graph representing a traffic network
//            with one source and one sink, we have for each arc an
//            associated capacity, base travel time and a
//            sensitivity. Travel time along a specific arc increases
//            as the flow approaches the capacity.
//
//            Given a fixed inflow we now wish to find the
//            configuration that minimizes the average travel time.


using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class TrafficNetworkModel : Model
  {

    public static double[] solve
    ( int      numberOfNodes,
      int      source_idx,
      int      sink_idx,
      int[]    arc_i,
      int[]    arc_j,
      double[] arcSensitivity,
      double[] arcCapacity,
      double[] arcBaseTravelTime,
      double   T)

    {
      using (Model M = new Model("Traffic Network"))
      {
        int n = numberOfNodes;        
        int m = arc_j.Length;

        double[] n_ones = new double[m]; for (int i = 0; i < m; ++i) n_ones[i] = 1.0;
        var NxN = Set.Make(n, n);
        Matrix basetime =
          Matrix.Sparse(n, n, arc_i, arc_j, arcBaseTravelTime);
        Matrix e =
          Matrix.Sparse(n, n, arc_i, arc_j, n_ones);
        Matrix e_e =
          Matrix.Sparse(n, n,
                        new int[]  { sink_idx}, new int[] { source_idx},
                        new double[] { 1.0 });

        int[,] sparsity = new int[m+1,2]; for (int i = 0; i < m; ++i) { sparsity[i,0] = arc_i[i]; sparsity[i,1] = arc_j[i]; }
        sparsity[m,0] = sink_idx;
        sparsity[m,1] = source_idx;

        Console.WriteLine("sparsity dims : {0} x {1}",sparsity.GetLength(0), sparsity.GetLength(1));
        
        double[] cs_inv = new double[m];
        double[] s_inv  = new double[m];
        for (int i = 0; i < m; ++i) cs_inv[i] = 1.0 / (arcSensitivity[i] * arcCapacity[i]);
        for (int i = 0; i < m; ++i) s_inv[i]  = 1.0 / arcSensitivity[i];
        Matrix cs_inv_matrix = Matrix.Sparse(n, n, arc_i, arc_j, cs_inv);
        Matrix s_inv_matrix  = Matrix.Sparse(n, n, arc_i, arc_j, s_inv);

        Variable x = M.Variable("traffic_flow", NxN, Domain.GreaterThan(0.0).Sparse(sparsity));
        Variable d = M.Variable("d",            NxN, Domain.GreaterThan(0.0).Sparse(sparsity));
        Variable z = M.Variable("z",            NxN, Domain.GreaterThan(0.0).Sparse(sparsity));


        // Set the objective:       
        M.Objective("Average travel time", ObjectiveSense.Minimize,
                    Expr.Mul(1.0 / T, Expr.Add(Expr.Dot(basetime, x), Expr.Sum(d))));

        // Set up constraints
        // Constraint (1a)
        M.Constraint("(1a)", Expr.Hstack(d.Pick(sparsity),
                                         z.Pick(sparsity),
                                         x.Pick(sparsity)).Slice(new int[]{0,0},new int[]{m,3} ), Domain.InRotatedQCone());
        
        // Constraint (1b)
        M.Constraint("(1b)",
                     Expr.Sub(Expr.Add(Expr.MulElm(z, e),
                                       Expr.MulElm(x, cs_inv_matrix)),
                              s_inv_matrix),
                     Domain.EqualsTo(0.0));

        // Constraint (2)
        M.Constraint("(2)", Expr.Sub(Expr.Add(Expr.MulDiag(x, e.Transpose()),
                                              Expr.MulDiag(x, e_e.Transpose())),
                                     Expr.Add(Expr.MulDiag(x.Transpose(), e),
                                              Expr.MulDiag(x.Transpose(), e_e))),
                     Domain.EqualsTo(0.0));
        // Constraint (3)
        M.Constraint("(3)", x.Index(sink_idx, source_idx), Domain.EqualsTo(T));

        M.Solve();
        M.WriteTask("TrafficNetworkModel.task");

        return x.Level();
      }
    }

    public static void Main(string[] args)
    {
      int n             = 4;
      int[]    arc_i    = new int[]    {  0,    0,    2,    1,    2 };
      int[]    arc_j    = new int[]    {  1,    2,    1,    3,    3 };
      double[] arc_base = new double[] {  4.0,  1.0,  2.0,  1.0,  6.0 };
      double[] arc_cap  = new double[] { 10.0, 12.0, 20.0, 15.0, 10.0 };
      double[] arc_sens = new double[] {  0.1,  0.7,  0.9,  0.5,  0.1 };

      double   T          = 20.0;
      int      source_idx = 0;
      int      sink_idx   = 3;

      try {
          var flow = solve(n, source_idx, sink_idx,
                           arc_i, arc_j,
                           arc_sens,
                           arc_cap,
                           arc_base,
                           T);
          System.Console.WriteLine("Optimal flow: {0} {1}",arc_i.Length,flow.Length);
          
          for (int i = 0; i < arc_i.Length; ++i)
            System.Console.WriteLine("\tflow node {0}->node {1} = {2}",
                                     arc_i[i], arc_j[i], flow[arc_i[i] * n + arc_j[i]]);
      }
      catch (mosek.fusion.SolutionError e )
        {
          Console.WriteLine("SolutionError: {0}",e.ToString());
          Console.WriteLine("{0}",e.StackTrace);
          throw;
        }
    }
  }
}

alan.cs

Listing 16.2 alan.cs Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:     alan.cs
//
//  Purpose: This file contains an implementation of the alan.gms (as
//  found in the GAMS online model collection) using Fusion.
//
//  The model is a simple portfolio choice model. The objective is to
//  invest in a number of assets such that we minimize the risk, while
//  requiring a certain expected return.
//
//  We operate with 4 assets (hardware,software, show-biz and treasure
//  bill). The risk is defined by the covariance matrix
//    Q = [[  4.0, 3.0, -1.0, 0.0 ],
//         [  3.0, 6.0,  1.0, 0.0 ],
//         [ -1.0, 1.0, 10.0, 0.0 ],
//         [  0.0, 0.0,  0.0, 0.0 ]]
//
//
//  We use the form Q = U^T * U, where U is a Cholesky factor of Q.
//

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class alan
  {
    /////////////////////////////////////////////////////////////////////
    // Problem data.

    // Security names
    private static string[]
    securities = { "hardware", "software", "show-biz", "t-bills" };
    // Two examples of mean returns on securities
    private static double[]
    mean1      =  {        9.0,        7.0,       11.0,       5.0 };
    private static double[]
    mean2       = {        8.0,        9.0,       12.0,       7.0 };
    // Target mean return
    private static double
    target     = 10.0;

    // Factor of covariance matrix.
    private static Matrix U =
      Matrix.Dense(
        new double[,]
    { {  2.0       ,  1.5       , -0.5       , 0.0 },
      {  0.0       ,  1.93649167,  0.90369611, 0.0 },
      {  0.0       ,  0.0       ,  2.98886824, 0.0 },
      {  0.0       ,  0.0       ,  0.0       , 0.0 }
    });
    private static int numsec = securities.Length;
    public static void Main(String[] args)
    {
      // Create a parametrized model
      using (Model M = new Model("alan"))
      {
        Variable x = M.Variable("x",        numsec, Domain.GreaterThan(0.0));
        Variable t = M.Variable("variance", Domain.GreaterThan(0.0));
        M.Objective("minvar", ObjectiveSense.Minimize, t.AsExpr());

        // sum securities to 1.0
        M.Constraint("wealth",  Expr.Sum(x), Domain.EqualsTo(1.0));
        // define target expected return
        Parameter mean = M.Parameter(numsec);
        M.Constraint("dmean", Expr.Dot(mean, x), Domain.GreaterThan(target));

        M.Constraint("q",Expr.Vstack(Expr.ConstTerm(1, 0.5),
                                 t.AsExpr(),
                                 Expr.Mul(U, x)),
                     Domain.InRotatedQCone());

        // Solve two instances of the problem
        solve(M, mean, x, mean1);
        solve(M, mean, x, mean2);
      }
    }

    public static void solve(Model M, Parameter mean, Variable x, double[] meanVal) {
      Console.Write("Solve with mean = ");
      for (int i = 0; i < numsec; ++i)
        Console.Write("  {0}", meanVal[i]);
      Console.WriteLine();

      mean.SetValue(meanVal);
      M.Solve();        

      double[] solx = x.Level();

      Console.Write("Solution = ");
      for (int i = 0; i < numsec; ++i)
        Console.Write(" {0}", solx[i]);
      Console.WriteLine("");        
    }
  }
}

baker.cs

Listing 16.3 baker.cs Click here to download.
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      baker.cs
//
// Purpose: Demonstrates a small linear problem.
//
// Source: "Linaer Algebra" by Knut Sydsaeter and Bernt Oeksendal.
//
// The problem: A baker has 150 kg flour, 22 kg sugar, 25 kg butter and two
// recipes:
//   1) Cakes, requiring 3.0 kg flour, 1.0 kg sugar and 1.2 kg butter per dozen.
//   2) Breads, requiring 5.0 kg flour, 0.5 kg sugar and 0.5 kg butter per dozen.
// Let the revenue per dozen cakes be $4 and the revenue per dozen breads be $6.
//
// We now wish to compute the combination of cakes and breads that will optimize
// the total revenue.

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class baker
  {
    private static string[] ingredientnames = { "Flour", "Sugar", "Butter" };
    private static double[] stock = { 150.0,   22.0,    25.0 };

    private static double[,] recipe_data =
    { { 3.0, 5.0 },
      { 1.0, 0.5 },
      { 1.2, 0.5 }
    };
    private static string[] productnames = { "Cakes", "Breads" };

    private static double[] revenue = { 4.0, 6.0 };

    public static void Main(string[] args)
    {
      Matrix recipe = Matrix.Dense(recipe_data);
      using (Model M = new Model("Recipe"))
      {
        // "production" defines the amount of each product to bake.
        Variable production = M.Variable("production",
                                         Set.Make(productnames),
                                         Domain.GreaterThan(0.0));
        // The objective is to maximize the total revenue.
        M.Objective("revenue",
                    ObjectiveSense.Maximize,
                    Expr.Dot(revenue, production));

        // The prodoction is constrained by stock:
        M.Constraint(Expr.Mul(recipe, production), Domain.LessThan(stock));
        M.SetLogHandler(Console.Out);

        // We solve and fetch the solution:
        M.Solve();
        double[] res = production.Level();
        Console.WriteLine("Solution:");
        for (int i = 0; i < res.Length; ++i)
        {
          Console.WriteLine(" Number of {0} : {1}", productnames[i], res[i]);
        }
        Console.WriteLine(" Revenue : ${0}",
                          res[0] * revenue[0] + res[1] * revenue[1]);
      }
    }
  }
}

breaksolver.cs

Listing 16.4 breaksolver.cs Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      breaksolver.cs
//
//  Purpose: Show how to break a long-running task.
//
//  Requires a parameter defining a timeout in seconds.
//


using System;
using System.Threading;
using mosek.fusion;
using System.Diagnostics;

namespace mosek.fusion.example
{
  public class breaksolver
  {
    private static void nshuffle(Random R, int[] a, int n)
    {
      for (int i = 0; i < n; ++i)
      {
        int idx = R.Next(i, a.Length);
        int tmp = a[i]; a[i] = a[idx]; a[idx] = tmp;
      }
    }

    public static void Main(string[] args)
    {
      long timeout = 5;

      int n = 200;    // number of binary variables
      int m = n / 5; // number of constraints
      int p = n / 5; // Each constraint picks p variables and requires that half of them are 1
      Random R = new Random(1234);
      using (Model M = new Model("SolveBinary"))
      {
        M.SetLogHandler(System.Console.Out);

        //Variable x = M.Variable("x", n, Domain.InRange(0,1));
        Variable x = M.Variable("x", n, Domain.Binary());
        M.Objective(ObjectiveSense.Minimize, Expr.Sum(x));
        //M.SetSolverParam("numThreads",1);

        int[] idxs  = new int[n]; for (int i = 0; i < n; ++i) idxs[i] = i;
        int[] cidxs = new int[p];

        for (var i = 0; i < m; ++i)
        {
          nshuffle(R, idxs, p);
          Array.Copy(idxs, cidxs, p);
          M.Constraint(Expr.Sum(x.Pick(cidxs)), Domain.EqualsTo(p / 2));
        }

        var T = new Thread(new ThreadStart(M.Solve));
        T.Start();

        Stopwatch w = new Stopwatch(); w.Start();
        while (true)
        {
          if (w.ElapsedMilliseconds > timeout * 1000)
          {
            Console.WriteLine("Solver terminated due to timeout!");
            M.BreakSolver();
            T.Join();
            break;
          }
          if (! T.IsAlive) {
            Console.WriteLine("Solver terminated before anything happened!");
            T.Join();
            break;
          }
        }

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

callback.cs

Listing 16.5 callback.cs Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      callback.cs
//
//   Purpose:   To demonstrate how to use the progress
//              callback.
//
//              Use this script as follows:
//              callback psim
//              callback dsim
//              callback intpnt
//
//              The first argument tells which optimizer to use
//              i.e. psim is primal simplex, dsim is dual simplex
//              and intpnt is interior-point.
//

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  class myCallback : mosek.DataCallback
  {
    double maxtime;
    Model M;

    public myCallback( double maxtime_,
                       Model M_)
    {
      maxtime = maxtime_;
      M = M_;
    }

    public override int callback( callbackcode caller,
                                  double[]     douinf,
                                  int[]        intinf,
                                  long[]       lintinf )
    {
      double opttime = 0.0;
      int itrn;
      double pobj, dobj, stime;

      switch (caller)
      {
        case callbackcode.begin_intpnt:
          Console.WriteLine("Starting interior-point optimizer");
          break;
        case callbackcode.intpnt:
          itrn    = intinf[(int) iinfitem.intpnt_iter      ];
          pobj    = douinf[(int) dinfitem.intpnt_primal_obj];
          dobj    = douinf[(int) dinfitem.intpnt_dual_obj  ];
          stime   = douinf[(int) dinfitem.intpnt_time      ];
          opttime = douinf[(int) dinfitem.optimizer_time   ];

          Console.WriteLine("Iterations: {0,-3}", itrn);
          Console.WriteLine("  Elapsed: Time: {0,6:F2}({1:F2})", opttime, stime);
          Console.WriteLine("  Primal obj.: {0,-18:E6}  Dual obj.: {1,018:E6}e", pobj, dobj);
          break;
        case callbackcode.end_intpnt:
          Console.WriteLine("Interior-point optimizer finished.");
          break;
        case callbackcode.begin_primal_simplex:
          Console.WriteLine("Primal simplex optimizer started.");
          break;
        case callbackcode.update_primal_simplex:
          itrn    = intinf[(int) iinfitem.sim_primal_iter  ];
          pobj    = douinf[(int) dinfitem.sim_obj          ];
          stime   = douinf[(int) dinfitem.sim_time         ];
          opttime = douinf[(int) dinfitem.optimizer_time   ];

          Console.WriteLine("Iterations: {0,-3}}", itrn);
          Console.WriteLine("  Elapsed time: {0,6:F2}({1:F2})", opttime, stime);
          Console.WriteLine("  Obj.: {0,-18:E6}", pobj );
          break;
        case callbackcode.end_primal_simplex:
          Console.WriteLine("Primal simplex optimizer finished.");
          break;
        case callbackcode.begin_dual_simplex:
          Console.WriteLine("Dual simplex optimizer started.");
          break;
        case callbackcode.update_dual_simplex:
          itrn    = intinf[(int) iinfitem.sim_dual_iter    ];
          pobj    = douinf[(int) dinfitem.sim_obj          ];
          stime   = douinf[(int) dinfitem.sim_time         ];
          opttime = douinf[(int) dinfitem.optimizer_time   ];
          Console.WriteLine("Iterations: {0,-3}}", itrn);
          Console.WriteLine("  Elapsed time: {0,6:F2}({1:F2})", opttime, stime);
          Console.WriteLine("  Obj.: {0,-18:E6}", pobj );
          break;
        case callbackcode.end_dual_simplex:
          Console.WriteLine("Dual simplex optimizer finished.");
          break;
        case callbackcode.begin_bi:
          Console.WriteLine("Basis identification started.");
          break;
        case  callbackcode.end_bi:
          Console.WriteLine("Basis identification finished.");
          break;
        default:
          break;
      }
      if (opttime >= maxtime)
      {
        Console.WriteLine("B too much time. Terminate it.");
        return 1;
      }
      return 0;
    }
  }

  class myProgressCallback : mosek.Progress
  {
    public override int progressCB( callbackcode caller )
    {
      //Handle caller code here
      return 0;
    }
  }

  public class callback
  {
    public static void Main(string[] args)
    {
      string slvr     = "intpnt";

      if (args.Length < 1)
      {
        Console.WriteLine("Usage: callback ( psim | dsim | intpnt )");
      }

      if (args.Length >= 1) slvr = args[0];

      // We create a large linear problem
      int n = 150;
      int m = 700;
      double[] A = new double[m * n];
      double[] b = new double[m];
      double[] c = new double[n];

      Random rnd = new Random();
      for (int i = 0; i < m * n; i++) A[i] = rnd.NextDouble();
      for (int i = 0; i < m; i++)   b[i] = rnd.NextDouble();
      for (int i = 0; i < n; i++)   c[i] = rnd.NextDouble();

      double maxtime = 0.07;

      Model M = new Model("callback");

      Variable x = M.Variable(n, Domain.Unbounded());
      M.Constraint(Expr.Mul(Matrix.Dense(m, n, A), x), Domain.LessThan(b));
      M.Objective(ObjectiveSense.Maximize, Expr.Dot(c, x));

      if ( slvr == "psim")
        M.SetSolverParam("optimizer", "primalSimplex");
      else if ( slvr == "dsim")
        M.SetSolverParam("optimizer", "dualSimplex");
      else if ( slvr == "intpnt")
        M.SetSolverParam("optimizer", "intpnt");

      M.SetDataCallbackHandler( new myCallback(maxtime, M) );

      M.Solve();
    }
  }
}

ceo1.cs

Listing 16.6 ceo1.cs Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      ceo1.cs
//
//   Purpose: Demonstrates how to solve the problem
//
//   minimize x1 + x2
//   such that
//            x1 + x2 + x3  = 1.0
//                x1,x2    >= 0.0
//   and      x1 >= x2 * exp(x3/x2)
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class ceo1
  {
    public static void Main(string[] args)
    {
      using (Model M = new Model("ceo1"))
      {

        Variable x = M.Variable("x", 3, Domain.Unbounded());

        // Create the constraint
        //      x[0] + x[1] + x[2] = 1.0
        M.Constraint("lc", Expr.Sum(x), Domain.EqualsTo(1.0));

        // Create the exponential conic constraint
        Constraint expc = M.Constraint("expc", x, Domain.InPExpCone());

        // Set the objective function to (x[0] + x[1])
        M.Objective("obj", ObjectiveSense.Minimize, Expr.Sum(x.Slice(0,2)));

        // Solve the problem
        M.Solve();

        // Get the linear solution values
        double[] solx = x.Level();
        Console.WriteLine("x1,x2,x3 = {0}, {1}, {2}", solx[0], solx[1], solx[2]);

        // Get conic solution of expc
        double[] expclvl = expc.Level();
        double[] expcsn  = expc.Dual();
        
        Console.Write("expc levels = {0}", expclvl[0]);
        for (int i = 1; i < expclvl.Length; ++i)
          Console.Write(", {0}", expclvl[i]);
        Console.WriteLine();

        Console.Write("expc dual conic var levels = {0}", expcsn[0]);
        for (int i = 1; i < expcsn.Length; ++i)
          Console.Write(", {0}", expcsn[i]);
        Console.WriteLine();
      }
    }
  }
}

cqo1.cs

Listing 16.7 cqo1.cs Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      cqo1.cs
//
//   Purpose: Demonstrates how to solve the problem
//
//   minimize y1 + y2 + y3
//   such that
//            x1 + x2 + 2.0 x3  = 1.0
//                    x1,x2,x3 >= 0.0
//   and
//            (y1,x1,x2) in C_3,
//            (y2,y3,x3) in K_3
//
//   where C_3 and K_3 are respectively the quadratic and
//   rotated quadratic cone of size 3 defined as
//       C_3 = { z1,z2,z3 :      z1 >= sqrt(z2^2 + z3^2) }
//       K_3 = { z1,z2,z3 : 2 z1 z2 >= z3^2              }
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class cqo1
  {
    public static void Main(string[] args)
    {
      using (Model M = new Model("cqo1"))
      {

        Variable x = M.Variable("x", 3, Domain.GreaterThan(0.0));
        Variable y = M.Variable("y", 3, Domain.Unbounded());

        // Create the aliases
        //      z1 = [ y[0],x[0],x[1] ]
        //  and z2 = [ y[1],y[2],x[2] ]
        Variable z1 = Var.Vstack(y.Index(0),  x.Slice(0, 2));
        Variable z2 = Var.Vstack(y.Slice(1, 3), x.Index(2));

        // Create the constraint
        //      x[0] + x[1] + 2.0 x[2] = 1.0
        double[] aval = new double[] {1.0, 1.0, 2.0};
        M.Constraint("lc", Expr.Dot(aval, x), Domain.EqualsTo(1.0));

        // Create the constraints
        //      z1 belongs to C_3
        //      z2 belongs to K_3
        // where C_3 and K_3 are respectively the quadratic and
        // rotated quadratic cone of size 3, i.e.
        //                 z1[0] >= sqrt(z1[1]^2 + z1[2]^2)
        //  and  2.0 z2[0] z2[1] >= z2[2]^2
        Constraint qc1 = M.Constraint("qc1", z1.AsExpr(), Domain.InQCone());
        Constraint qc2 = M.Constraint("qc2", z2.AsExpr(), Domain.InRotatedQCone());

        // Set the objective function to (y[0] + y[1] + y[2])
        M.Objective("obj", ObjectiveSense.Minimize, Expr.Sum(y));

        // Solve the problem
        M.Solve();

        // Get the linear solution values
        double[] solx = x.Level();
        double[] soly = y.Level();
        Console.WriteLine("x1,x2,x3 = {0}, {1}, {2}", solx[0], solx[1], solx[2]);
        Console.WriteLine("y1,y2,y3 = {0}, {1}, {2}", soly[0], soly[1], soly[2]);

        // Get conic solution of qc1
        double[] qc1lvl = qc1.Level();
        double[] qc1sn  = qc1.Dual();
        
        Console.Write("qc1 levels = {0}", qc1lvl[0]);
        for (int i = 1; i < qc1lvl.Length; ++i)
          Console.Write(", {0}", qc1lvl[i]);
        Console.WriteLine();

        Console.Write("qc1 dual conic var levels = {0}", qc1sn[0]);
        for (int i = 1; i < qc1sn.Length; ++i)
          Console.Write(", {0}", qc1sn[i]);
        Console.WriteLine();
      }
    }
  }
}

diet.cs

Listing 16.8 diet.cs Click here to download.
/*
  File : diet.cs

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic linear model.
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class diet : Model
  {
    public class DataException : FusionException
    {
      public DataException(string msg) : base(msg) {  }
    }

    public double[] dailyAllowance;
    public Matrix nutrientValue;

    private Variable   dailyPurchase;
    private Constraint dailyNutrients;

    public diet(string    name,
                string[]  foods,
                string[]  nutrients,
                double[]  daily_allowance,
                double[,] nutritive_value)
    : base(name)
    {
      bool finished = false;
      try
      {
        dailyAllowance = daily_allowance;
        nutrientValue = (Matrix.Dense(nutritive_value)).Transpose();

        int M = foods.Length;
        int N = nutrients.Length;

        if (dailyAllowance.Length != N)
          throw new DataException(
            "Length of daily_allowance does not match the number of nutrients");
        if (nutrientValue.NumColumns() != M)
          throw new DataException(
            "Number of rows in nutrient_value does not match the number of foods");
        if (nutrientValue.NumRows() != N)
          throw new DataException(
            "Number of columns in nutrient_value does not match the number of nutrients");

        dailyPurchase = Variable("Daily Purchase",
                                 Set.Make(foods),
                                 Domain.GreaterThan(0.0).WithNamesOnAxis(foods,0));
        dailyNutrients =
          Constraint("Nutrient Balance",
                     Expr.Mul(nutrientValue, dailyPurchase),
                     Domain.GreaterThan(dailyAllowance).WithNamesOnAxis(nutrients,0));
        Objective(ObjectiveSense.Minimize, Expr.Sum(dailyPurchase));
        WriteTask("diet.ptf");
        finished = true;
      }
      finally
      {
        if (! finished)
          Dispose();
      }
    }

    public double[] getDailyPurchase()
    {
      return dailyPurchase.Level();
    }

    public double[] getDailyNutrients()
    {
      return dailyNutrients.Level();
    }


    /* Main class with data definitions */
    public static void Main(string[] argv)
    {
      string[] nutrient_unit = {
        "thousands",  "grams",        "grams",
        "milligrams", "thousand ius", "milligrams",
        "milligrams", "milligrams",   "milligrams"
      };
      string[] nutrients = {
        "calorie",    "protein",      "calcium",
        "iron",       "vitamin a",    "vitamin b1",
        "vitamin b2", "niacin",       "vitamin c"
      };
      string[] foods = {
        "wheat",        "cornmeal", "cannedmilk", "margarine", "cheese",
        "peanut butter", "lard",     "liver",      "porkroast", "salmon",
        "greenbeans",   "cabbage",  "onions",     "potatoes",  "spinach",
        "sweet potatos", "peaches",  "prunes",     "limabeans", "navybeans"
      };
      double[,] nutritive_value = {
        //  calorie       calcium      vitamin a        vitamin b2      vitamin c
        //         protein        iron           vitamin b1      niacin
        {44.7,  1411,   2.0,   365,    0,      55.4,   33.3,  441,     0},  // wheat
        {36,     897,   1.7,    99,   30.9,    17.4,    7.9,  106,     0},  // cornmeal
        { 8.4,   422,  15.1,     9,   26,       3,     23.5,   11,    60},  // cannedmilk
        {20.6,    17,   0.6,     6,   55.8,     0.2,    0,      0,     0},  // margarine
        { 7.4,   448,  16.4,    19,   28.1,     0.8,   10.3,    4,     0},  // cheese
        {15.7,   661,   1,      48,     0,      9.6,    8.1,  471,     0},  // peanut butter
        {41.7,     0,   0,       0,    0.2,     0,      0.5,    5,     0},  // lard
        { 2.2,   333,   0.2,   139,  169.2,     6.4,   50.8,  316,   525},  // liver
        { 4.4,   249,   0.3,    37,    0,      18.2,    3.6,   79,     0},  // porkroast
        { 5.8,   705,   6.8,    45,    3.5,     1,      4.9,  209,     0},  // salmon
        { 2.4,   138,   3.7,    80,   69,       4.3,    5.8,   37,   862},  // greenbeans
        { 2.6,   125,   4,      36,    7.2,     9,      4.5,   26,  5369},  // cabbage
        { 5.8,   166,   3.8,    59,   16.6,     4.7,    5.9,   21,  1184},  // onions
        {14.3,   336,   1.8,   118,    6.7,    29.4,    7.1,  198,  2522},  // potatoes
        { 1.1,   106,   0.0,   138,  918.4,     5.7,   13.8,   33,  2755},  // spinach
        { 9.6,   138,   2.7,    54,  290.7,     8.4,    5.4,   83,  1912},  // sweet potatos
        { 8.5,    87,   1.7,   173,   86.8,     1.2,    4.3,   55,    57},  // peaches
        {12.8,    99,   2.5,   154,   85.7,     3.9,    4.3,   65,   257},  // prunes
        {17.4,  1055,   3.7,   459,    5.1,    26.9,   38.2,   93,     0},  // limabeans
        {26.9,  1691,  11.4,   792,    0,      38.4,   24.6,  217,     0}
      };// navybeans
      double[] daily_allowance =
      {   3.0,    70.0, 0.8,    12.0,  5.0,     1.8,    2.7,   18.0,  75.0 };

      using (diet M = new diet("Stinglers Diet Problem",
                               foods,
                               nutrients,
                               daily_allowance,
                               nutritive_value))
      {
        M.Solve();

        try
        {
          Console.WriteLine("Solution (x 365):");
          double[] x = M.getDailyPurchase();
          for (int i = 0; i < x.Length; ++i)
            Console.WriteLine("  {0,15} : {1,10:0.00} $", foods[i], 365 * x[i]);

          Console.WriteLine("Nutrients (daily):");
          double[] y = M.getDailyNutrients();
          for (int i = 0; i < y.Length; ++i)
            Console.WriteLine("  {0,15} : {1,10:0.00} {2,-12} ({3})",
                              nutrients[i],
                              y[i],
                              nutrient_unit[i],
                              daily_allowance[i]);
        }
        catch (SolutionError e)
        {
          Console.WriteLine("Solution error : {0}", e.ToString());
        }
      }
    }
  }
}

djc1.cs

Listing 16.9 djc1.cs Click here to download.
////
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      djc1.cs
//
//   Purpose: Demonstrates how to solve the problem with two disjunctions:
//
//      minimize    2x0 + x1 + 3x2 + x3
//      subject to   x0 + x1 + x2 + x3 >= -10
//                  (x0-2x1<=-1 and x2=x3=0) or (x2-3x3<=-2 and x1=x2=0)
//                  x0=2.5 or x1=2.5 or x2=2.5 or x3=2.5
////
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class djc1
  {
    public static void Main(string[] args)
    {
      Model M = new Model("djc1");

      // Create variable 'x' of length 4
      Variable x = M.Variable("x", 4);

      // First disjunctive constraint
      M.Disjunction( DJC.AND( DJC.Term(Expr.Dot(new double[]{1,-2,0,0}, x), Domain.LessThan(-1)), // x0 - 2x1 <= -1  
                              DJC.Term(x.Slice(2, 4), Domain.EqualsTo(0)) ),                      // x2 = x3 = 0
                     DJC.AND( DJC.Term(Expr.Dot(new double[]{0,0,1,-3}, x), Domain.LessThan(-2)), // x2 - 3x3 <= -2
                              DJC.Term(x.Slice(0, 2), Domain.EqualsTo(0)) ) );                    // x0 = x1 = 0

      // Second disjunctive constraint
      // Array of terms reading x_i = 2.5 for i = 0,1,2,3
      Term[] terms = new Term[4];
      for(int i = 0; i < 4; i++)
        terms[i] = DJC.Term(x.Index(i), Domain.EqualsTo(2.5));
      // The disjunctive constraint from the array of terms
      M.Disjunction(terms);

      // The linear constraint
      M.Constraint(Expr.Sum(x), Domain.GreaterThan(-10));

      // Objective
      M.Objective(ObjectiveSense.Minimize, Expr.Dot(new double[]{2,1,3,1}, x));

      // Useful for debugging
      M.WriteTask("djc.ptf");         // Save to a readable file
      M.SetLogHandler(Console.Out);   // Enable log output

      // Solve the problem
      M.Solve();

      // Get the solution values
      if (M.GetPrimalSolutionStatus() == SolutionStatus.Optimal) {
        double[] sol = x.Level();
        Console.WriteLine("[x0,x1,x2,x3] = [{0}, {1}, {2}, {3} ]", sol[0], sol[1], sol[2], sol[3]);
      }
      else {
        Console.WriteLine("Another solution status");
      }
    }
  }
}

duality.cs

Listing 16.10 duality.cs Click here to download.
/*
  File : duality.cs

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Shows how to read the dual value of a constraint
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class duality
  {
    public static void Main(string[] args)
    {
      double[,] A =  new double[,] { { -0.5, 1.0 }  };
      double[]  b =  new double[] { 1.0 };
      double[]  c =  new double[] { 1.0, 1.0 };

      using (Model M = new Model("duality"))
      {
        Variable x = M.Variable("x", 2, Domain.GreaterThan(0.0));

        Constraint con = M.Constraint(Expr.Sub(b, Expr.Mul(Matrix.Dense(A), x)), Domain.EqualsTo(0.0));

        M.Objective("obj", ObjectiveSense.Minimize, Expr.Dot(c, x));

        M.Solve();
        double[] xsol = x.Level();
        double[] ysol = con.Dual();

        Console.WriteLine("x1,x2,y = %{0}, %{1}, %{2}\n", xsol[0], xsol[1], ysol[0]);
      }
    }
  }
}

elastic.cs

Listing 16.11 elastic.cs Click here to download.
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      elastic.cs
//
// Purpose: Demonstrates model parametrization on the example of an elastic net linear regression:
//
//          min_x  |Ax-b|_2 + lambda1*|x|_1 + lambda2*|x|_2

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class elastic
  {
    // Construct the model with parameters b, lambda1, lambda2
    // and with prescribed matrix A
    public static Model initializeModel(int m, int n, double[,] A) {
      Model M = new Model();
      Variable x = M.Variable("x", n);

      // t >= |Ax-b|_2 where b is a parameter
      Parameter b = M.Parameter("b", m);
      Variable  t = M.Variable();
      M.Constraint(Expr.Vstack(t, Expr.Sub(Expr.Mul(A, x), b)), Domain.InQCone());

      // p_i >= |x_i|, i=1..n
      Variable p = M.Variable(n);
      M.Constraint(Expr.Hstack(p, x), Domain.InQCone());

      // q >= |x|_2
      Variable q = M.Variable();
      M.Constraint(Expr.Vstack(q, x), Domain.InQCone());

      // Objective, parametrized with lambda1, lambda2
      // t + lambda1*sum(p) + lambda2*q
      Parameter lambda1 = M.Parameter("lambda1");
      Parameter lambda2 = M.Parameter("lambda2");
      Expression obj = Expr.Add(new Expression[] {t, Expr.Mul(lambda1, Expr.Sum(p)), Expr.Mul(lambda2, q)});
      M.Objective(ObjectiveSense.Minimize, obj);

      // Return the ready model
      return M;
    }

    public static void smallExample() {
      //Create a small example
      int m = 4;
      int n = 2;
      double[,] A = { {1.0,   2.0},
                      {3.0,   4.0},
                      {-2.0, -1.0},
                      {-4.0, -3.0} };
      double[] sol;
      Model M = initializeModel(m, n, A);

      // For convenience retrieve some elements of the model
      Parameter b = M.GetParameter("b");
      Parameter lambda1 = M.GetParameter("lambda1");
      Parameter lambda2 = M.GetParameter("lambda2");
      Variable x = M.GetVariable("x");

      // First solve
      b.SetValue(new double[]{0.1, 1.2, -1.1, 3.0});
      lambda1.SetValue(0.1);
      lambda2.SetValue(0.01);

      M.Solve();
      sol = x.Level();
      Console.WriteLine("Objective {0}, solution {1}, {2}", M.PrimalObjValue(), sol[0], sol[1]);

      // Increase lambda1
      lambda1.SetValue(0.5);
      
      M.Solve();
      sol = x.Level();
      Console.WriteLine("Objective {0}, solution {1}, {2}", M.PrimalObjValue(), sol[0], sol[1]);

      // Now change the data completely
      b.SetValue(new double[] {1.0, 1.0, 1.0, 1.0});
      lambda1.SetValue(0.0);
      lambda2.SetValue(0.0);
      
      M.Solve();
      sol = x.Level();
      Console.WriteLine("Objective {0}, solution {1}, {2}", M.PrimalObjValue(), sol[0], sol[1]);

      // And increase lamda2
      lambda2.SetValue(1.4145);
      
      M.Solve();
      sol = x.Level();
      Console.WriteLine("Objective {0}, solution {1}, {2}", M.PrimalObjValue(), sol[0], sol[1]);

      M.Dispose();
    }

    public static void Main() {
      smallExample();
    }
  }
}

facility_location.cs

Listing 16.12 facility_location.cs Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      facility_location.cs
//
//  Purpose: Demonstrates a small one-facility location problem.
//
//  Given 10 customers placed in a grid we wish to place a facility
//  somewhere so that the total sum of distances to customers is
//  minimized.
//
//  The problem is formulated as a conic optimization problem as follows.
//  Let f=(fx,fy) be the (unknown) location of the facility, and let
//  c_i=(cx_i,cy_i) be the (known) customer locations; then we wish to
//  minimize
//      sum_i || f - c_i ||
//  where
//      ||.||
//  denotes the euclidian norm.
//  This is formulated as
//
//  minimize   sum(d_i)
//  such that  d_i ^ 2 > tx_i ^ 2 + ty_i ^ 2, for all i
//             tx_i = cx_i - fx, for all i
//             ty_i = cy_i - fy, for all i
//             d_i > 0, for all i
////

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class facility_location
  {
    // Customer locations
    private static Matrix  customerloc =
      Matrix.Dense (new double[,]
    { { 12,  2 },
      { 15, 13 },
      { 10,  8 },
      {  0, 10 },
      {  6, 13 },
      {  5,  8 },
      { 10, 12 },
      {  4,  6 },
      {  5,  2 },
      {  1, 10 }
    } );
    private static int N = customerloc.NumRows();
    public static void Main(string[] args)
    {
      using (Model M = new Model("FacilityLocation"))
      {
        // Variable holding the facility location
        Variable f = M.Variable("facility", Set.Make(1, 2), Domain.Unbounded());
        // Variable defining the euclidian distances to each customer
        Variable d = M.Variable("dist", Set.Make(N, 1), Domain.GreaterThan(0.0));
        // Variable defining the x and y differences to each customer;
        Variable t = M.Variable("t", Set.Make(N, 2), Domain.Unbounded());
        M.Constraint("dist measure",
        Var.Hstack(new Variable[] { d, t }),
        Domain.InQCone(N, 3));

        Variable fxy = Var.Repeat(f, N);
        M.Constraint("xy diff", Expr.Add(t, fxy), Domain.EqualsTo(customerloc));

        M.Objective("total_dist", ObjectiveSense.Minimize, Expr.Sum(d));

        M.Solve();
        M.WriteTask("facility_location.task");
        double[] floc = f.Level();
        Console.WriteLine("Facility location = {0},{1}", floc[0], floc[1]);
      }
    }
  }
}

gp1.cs

Listing 16.13 gp1.cs Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      gp1.cs
//
//   Purpose:   Demonstrates how to solve a simple Geometric Program (GP)
//              cast into conic form with exponential cones and log-sum-exp.
//
//              Example from
//                https://gpkit.readthedocs.io/en/latest/examples.html//maximizing-the-volume-of-a-box
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class gp1
  {
    // Models log(sum(exp(Ax+b))) <= 0.
    // Each row of [A b] describes one of the exp-terms
    public static void logsumexp(Model      M, 
                                 double[,]  A, 
                                 Variable   x,
                                 double[]   b)
    {
      int k = A.GetLength(0);
      Variable u = M.Variable(k);
      M.Constraint(Expr.Sum(u), Domain.EqualsTo(1.0));
      M.Constraint(Expr.Hstack(u,
                               Expr.ConstTerm(k, 1.0),
                               Expr.Add(Expr.Mul(A, x), b)), Domain.InPExpCone());
    }

    // maximize     h*w*d
    // subjecto to  2*(h*w + h*d) <= Awall
    //              w*d <= Afloor
    //              alpha <= h/w <= beta
    //              gamma <= d/w <= delta
    //
    // Variable substitutions:  h = exp(x), w = exp(y), d = exp(z).
    //
    // maximize     x+y+z
    // subject      log( exp(x+y+log(2/Awall)) + exp(x+z+log(2/Awall)) ) <= 0
    //                              y+z <= log(Afloor)
    //              log( alpha ) <= x-y <= log( beta )
    //              log( gamma ) <= z-y <= log( delta )
    public static double[] max_volume_box(double Aw, double Af, 
                                          double alpha, double beta, double gamma, double delta)
    {
      using (Model M = new Model("max_vol_box"))
      {
        Variable xyz = M.Variable(3);
        M.Objective("Objective", ObjectiveSense.Maximize, Expr.Sum(xyz));
        
        logsumexp(M, 
                  new double[,] {{1,1,0}, {1,0,1}}, 
                  xyz, 
                  new double[] {System.Math.Log(2.0/Aw), System.Math.Log(2.0/Aw)});
        
        M.Constraint(Expr.Dot(new double[] {0,1,1}, xyz), Domain.LessThan(System.Math.Log(Af)));
        M.Constraint(Expr.Dot(new double[] {1,-1,0}, xyz), Domain.InRange(System.Math.Log(alpha),System.Math.Log(beta)));
        M.Constraint(Expr.Dot(new double[] {0,-1,1}, xyz), Domain.InRange(System.Math.Log(gamma),System.Math.Log(delta)));
        
        M.SetLogHandler(Console.Out);
        M.Solve();
        
        double[] xyzVal = xyz.Level();
        double[] hwdVal = new double[3];
        for(int i=0; i<3; i++) hwdVal[i] = System.Math.Exp(xyzVal[i]);

        return hwdVal;
      }
    }

    public static void Main(String[] args)
    {
      double Aw    = 200.0;
      double Af    = 50.0;
      double alpha = 2.0;
      double beta  = 10.0;
      double gamma = 2.0;
      double delta = 10.0;
      
      double[] hwd = max_volume_box(Aw, Af, alpha, beta, gamma, delta);

      Console.WriteLine("h={0:f4} w={1:f4} d={2:f4}", hwd[0], hwd[1], hwd[2]);
    }
  }
}

lo1.cs

Listing 16.14 lo1.cs Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      lo1.cs
//
//  Purpose: Demonstrates how to solve the problem
//
//  maximize 3*x0 + 1*x1 + 5*x2 + x3
//  such that
//           3*x0 + 1*x1 + 2*x2        = 30,
//           2*x0 + 1*x1 + 3*x2 + 1*x3 > 15,
//                  2*x1 +      + 3*x3 < 25
//  and
//           x0,x1,x2,x3 > 0,
//           0 < x1 < 10
////
using System;
using mosek.fusion;
namespace mosek.fusion.example
{
  public class lo1
  {
    public static void Main(string[] args)
    {
      double[][] A =
      { new double[] { 3.0, 1.0, 2.0, 0.0 },
        new double[] { 2.0, 1.0, 3.0, 1.0 },
        new double[] { 0.0, 2.0, 0.0, 3.0 }
      };
      double[] c = { 3.0, 1.0, 5.0, 1.0  };

      // Create a model with the name 'lo1'
      Model M = new Model("lo1");

      // Create variable 'x' of length 4
      Variable x = M.Variable("x", 4, Domain.GreaterThan(0.0));

      // Create constraints
      M.Constraint(x.Index(1), Domain.LessThan(10.0));
      
      M.Constraint("c1", Expr.Dot(A[0], x), Domain.EqualsTo(30.0));
      M.Constraint("c2", Expr.Dot(A[1], x), Domain.GreaterThan(15.0));
      M.Constraint("c3", Expr.Dot(A[2], x), Domain.LessThan(25.0));

      // Set the objective function to (c^t * x)
      M.Objective("obj", ObjectiveSense.Maximize, Expr.Dot(c, x));

      // Solve the problem
      M.Solve();

      // Get the solution values
      double[] sol = x.Level();
      Console.WriteLine("[x0,x1,x2,x3] = [{0}, {1}, {2}, {3} ]", sol[0], sol[1], sol[2], sol[3]);
    }
  }
}

logistic.cs

Listing 16.15 logistic.cs Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      logistic.cs
//
// Purpose: Implements logistic regression with regulatization.
//
//          Demonstrates using the exponential cone and log-sum-exp in Fusion.

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class logistic
  {
    // t >= log( 1 + exp(u) ) coordinatewise
    public static void softplus(Model      M,
                                Expression t,
                                Expression u)
    {
      int n = t.GetShape()[0];
      Variable z1 = M.Variable(n);
      Variable z2 = M.Variable(n);
      M.Constraint(Expr.Add(z1, z2), Domain.EqualsTo(1));
      M.Constraint(Expr.Hstack(z1, Expr.ConstTerm(n, 1.0), Expr.Sub(u,t)), Domain.InPExpCone());
      M.Constraint(Expr.Hstack(z2, Expr.ConstTerm(n, 1.0), Expr.Neg(t)), Domain.InPExpCone());
    }

    // Model logistic regression (regularized with full 2-norm of theta)
    // X - n x d matrix of data points
    // y - length n vector classifying training points
    // lamb - regularization parameter
    public static Model logisticRegression(double[,]  X, 
                                           bool[]     y,
                                           double     lamb)
    {
      int n = X.GetLength(0);
      int d = X.GetLength(1);       // num samples, dimension
      
      Model M = new Model();   

      Variable theta = M.Variable("theta", d);
      Variable t     = M.Variable(n);
      Variable reg   = M.Variable();

      M.Objective(ObjectiveSense.Minimize, Expr.Add(Expr.Sum(t), Expr.Mul(lamb,reg)));
      M.Constraint(Var.Vstack(reg, theta), Domain.InQCone());
      
      double[] signs = new double[n];
      for(int i = 0; i < n; i++)
        if (y[i]) signs[i] = -1; else signs[i] = 1;

      softplus(M, t, Expr.MulElm(Expr.Mul(X, theta), signs));

      return M;
    }

    public static void Main(String[] args)
    {
      // Test: detect and approximate a circle using degree 2 polynomials
      int n = 30;
      double[,] X = new double[n*n, 6];
      bool[] Y     = new bool[n*n];

      for(int i=0; i<n; i++) 
      for(int j=0; j<n; j++)
      {
        int k = i*n+j;
        double x = -1 + 2.0*i/(n-1);
        double y = -1 + 2.0*j/(n-1);
        X[k,0] = 1.0; X[k,1] = x; X[k,2] = y; X[k,3] = x*y;
        X[k,4] = x*x; X[k,5] = y*y;
        Y[k] = (x*x+y*y>=0.69);
      }

      Model M = logisticRegression(X, Y, 0.1);
      Variable theta = M.GetVariable("theta");

      M.Solve();
      for(int i=0; i<6; i++)
        Console.WriteLine(theta.Level()[i]);
    }
  }
}

lownerjohn_ellipsoid.cs

Listing 16.16 lownerjohn_ellipsoid.cs Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      lownerjohn_ellipsoid.cs

   Purpose:
   Computes the Lowner-John inner and outer ellipsoidal
   approximations of a polytope.


   The inner ellipsoidal approximation to a polytope

      S = { x \in R^n | Ax < b }.

   maximizes the volume of the inscribed ellipsoid,

      { x | x = C*u + d, || u ||_2 <= 1 }.

   The volume is proportional to det(C)^(1/n), so the
   problem can be solved as

     maximize         t
     subject to       t       <= det(C)^(1/n)
                 || C*ai ||_2 <= bi - ai^T * d,  i=1,...,m
                   C is PSD

   which is equivalent to a mixed conic quadratic and semidefinite
   programming problem.


   The outer ellipsoidal approximation to a polytope given
   as the convex hull of a set of points

       S = conv{ x1, x2, ... , xm }

   minimizes the volume of the enclosing ellipsoid,

     { x | || P*(x-c) ||_2 <= 1 }

   The volume is proportional to det(P)^{-1/n}, so the problem can
   be solved as

     minimize         t
     subject to       t       >= det(P)^(-1/n)
                 || P*xi + c ||_2 <= 1,  i=1,...,m
                   P is PSD.

   References:
   [1] "Lectures on Modern Optimization", Ben-Tal and Nemirovski, 2000.
 */

using System;

namespace mosek.fusion.example
{
  public class lownerjohn_ellipsoid
  {

    /**
        Purpose: Models the hypograph of the n-th power of the
        determinant of a positive definite matrix.

        The convex set (a hypograph)

          C = { (X, t) \in S^n_+ x R |  t <= det(X)^{1/n} },

        can be modeled as the intersection of a semidefinite cone

          [ X, Z; Z^T Diag(Z) ] >= 0

        and a geometric mean cone

          t <= (Z11*Z22*...*Znn)^{1/n}

        References:
        [1] "Lectures on Modern Optimization", Ben-Tal and Nemirovski, 2000.
     */
    public static Variable det_rootn(Model M, Variable t, int n)
    {
      // Setup variables
      Variable Y = M.Variable(Domain.InPSDCone(2 * n));

      Variable X   = Y.Slice(new int[]{0, 0}, new int[]{n, n});
      Variable Z   = Y.Slice(new int[]{0, n}, new int[]{n, 2 * n});
      Variable DZ  = Y.Slice(new int[]{n, n}, new int[]{2 * n, 2 * n});

      // Z is lower-triangular
      int[,] low_tri = new int[n*(n-1)/2,2];
      int k = 0;
      for(int i = 0; i < n; i++)
        for(int j = i+1; j < n; j++)
          { low_tri[k,0] = i; low_tri[k,1] = j; ++k; }
      M.Constraint(Z.Pick(low_tri), Domain.EqualsTo(0.0));
      // DZ = Diag(Z)
      M.Constraint(Expr.Sub(DZ, Expr.MulElm(Z, Matrix.Eye(n))), Domain.EqualsTo(0.0));

      // t^n <= (Z11*Z22*...*Znn)
      M.Constraint(Expr.Vstack(DZ.Diag(), t), Domain.InPGeoMeanCone());

      // Return an n x n PSD variable which satisfies t <= det(X)^(1/n)
      return X;
    }

    public static Tuple<double[], double[]> lownerjohn_inner(double[][] A, double[] b)
    {
      using( Model M = new Model("lownerjohn_inner"))
      {
        int m = A.Length;
        int n = A[0].Length;

        // Setup variables
        Variable t = M.Variable("t", 1, Domain.GreaterThan(0.0));
        Variable C = det_rootn(M, t, n);
        Variable d = M.Variable("d", n, Domain.Unbounded());

        // (bi - ai^T*d, C*ai) \in Q
        for (int i = 0; i < m; ++i)
          M.Constraint("qc" + i, Expr.Vstack(Expr.Sub(b[i], Expr.Dot(A[i], d)), Expr.Mul(C, A[i])),
                       Domain.InQCone().Axis(0) );

        // Objective: Maximize t
        M.Objective(ObjectiveSense.Maximize, t);
        M.Solve();

        return Tuple.Create(C.Level(), d.Level());
      }
    }

    public static Tuple<double[], double[]> lownerjohn_outer(double[,] x)
    {
      using( Model M = new Model("lownerjohn_outer") )
      {
        int m = x.GetLength(0);
        int n = x.GetLength(1);

        // Setup variables
        Variable t = M.Variable("t", 1, Domain.GreaterThan(0.0));
        Variable P = det_rootn(M, t, n);
        Variable c = M.Variable("c", Domain.Unbounded().WithShape(1,n));

        // (1, P(*xi+c)) \in Q
        M.Constraint("qc",
                     Expr.Hstack(Expr.Ones(m), Expr.Sub(Expr.Mul(x,P), Expr.Repeat(c,m,0))),
                     Domain.InQCone().Axis(1));

        // Objective: Maximize t
        M.Objective(ObjectiveSense.Maximize, t);
        M.Solve();

        return Tuple.Create(P.Level(), c.Level());
      }
    }

    public static void Main(string[] argv)
    {
      double[,] p = new double[,] { {0.0, 0.0},
                                    {1.0, 3.0},
                                    {5.5, 4.5},
                                    {7.0, 4.0},
                                    {7.0, 1.0},
                                    {3.0, -2.0} };
      int n = p.GetLength(0);
      double[][] A = new double[n][];
      double[]   b = new double[n];
      for (int i = 0; i < n; ++i)
      {
        A[i] = new double[] { -p[i,1] + p[(i - 1 + n) % n,1], p[i,0] - p[(i - 1 + n) % n,0] };
        b[i] = A[i][0] * p[i,0] + A[i][1] * p[i,1];
      }

      Tuple<double[], double[]> ResInn = lownerjohn_inner(A, b);
      Tuple<double[], double[]> ResOut = lownerjohn_outer(p);
      double[] Ci = ResInn.Item1;
      double[] di = ResInn.Item2;
      double[] Po = ResOut.Item1;
      double[] co = ResOut.Item2;
      Console.WriteLine("Inner:");
      Console.WriteLine("  C = [{0}]", string.Join(", ", Ci));
      Console.WriteLine("  d = [{0}]", string.Join(", ", di));
      Console.WriteLine("Outer:");
      Console.WriteLine("  P = [{0}]", string.Join(", ", Po));
      Console.WriteLine("  c = [{0}]", string.Join(", ", co));
    }
  }
}

lpt.cs

Listing 16.17 lpt.cs Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      lpt.cs
//
// Purpose:  Demonstrates how to solve the multi-processor
//           scheduling problem using the Fusion API.
//


using System;
using System.Collections.Generic;

using mosek.fusion;

namespace mosek.fusion.example
{
  public class lpt
  {
    public static void Main(string [] args)
    {
      int n = 30;           //Number of tasks
      int m = 6;              //Number of processors

      double lb = 1.0;         //The range of lengths of short tasks
      double ub = 5.0;

      double sh = 0.8;        //The proportion of short tasks
      int n_short = (int)(n * sh);
      int n_long = n - n_short;

      double[] T = new double[n];
      Random gen = new Random(0);
      for (int i = 0; i < n_short; i++)
        T[i] = gen.NextDouble() * (ub - lb) + lb;
      for (int i = n_short; i < n; i++)
        T[i] = 20 * (gen.NextDouble() * (ub - lb) + lb);
      Array.Sort<double>(T);

      Model M = new Model("Multi-processor scheduling");

      Variable x = M.Variable("x", new int[] {m, n}, Domain.Binary());
      Variable t = M.Variable("t", 1, Domain.Unbounded());

      M.Constraint( Expr.Sum(x, 0), Domain.EqualsTo(1.0) );
      M.Constraint( Expr.Sub( Var.Repeat(t, m), Expr.Mul(x, T)), Domain.GreaterThan(0.0) );

      M.Objective( ObjectiveSense.Minimize, t );

      //Computing LPT solution
      double [] init    = new double[n * m];
      double [] schedule = new double[m];

      for (int i = n - 1; i >= 0; i--)
      {
        int next = 0;
        for (int j = 1; j < m; j++)
          if (schedule[j] < schedule[next]) next = j;
        schedule[next] += T[i];
        init[next * n + i] = 1;
      }

      //Comment this line to switch off feeding in the initial LPT solution
      x.SetLevel(init);

      M.SetLogHandler(Console.Out);

      M.SetSolverParam("mioTolRelGap", 0.01);
      M.Solve();

      try
      {
        Console.Write("initial solution: \n");
        for (int i = 0; i < m; i++)
        {
          Console.Write("M {0} [", i);
          for (int y = 0; y < n; y++)
            Console.Write( "{0}, ", (int) init[i * n + y]);
          Console.Write("]\n");
        }
        Console.Write("MOSEK solution:\n");
        for (int i = 0; i < m; i++)
        {
          Console.Write("M {0} [", i);
          for (int y = 0; y < n; y++)
            Console.Write( "{0}, ", (int)x.Index(i, y).Level()[0]);
          Console.Write("]\n");
        }
      }
      catch (SolutionError e) {}
    }
  }
}

mico1.cs

Listing 16.18 mico1.cs Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    mico1.cs
//
//    Purpose:  Demonstrates how to solve a small mixed
//              integer conic optimization problem.
//
//              minimize    x^2 + y^2
//              subject to  x >= e^y + 3.8
//                          x, y - integer
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class mico1
  {
    public static void Main(string[] args)
    {
      using (Model M = new Model("mico1"))
      {
        Variable x = M.Variable(Domain.Integral(Domain.Unbounded()));
        Variable y = M.Variable(Domain.Integral(Domain.Unbounded()));
        Variable t = M.Variable();

        M.Constraint(Expr.Vstack(t, x, y), Domain.InQCone());
        M.Constraint(Expr.Vstack(Expr.Sub(x, 3.8), 1, y), Domain.InPExpCone());

        M.Objective(ObjectiveSense.Minimize, t);

        M.Solve();

        Console.WriteLine("x, y = {0}, {1}", x.Level()[0], y.Level()[0]);
      }
    }
  }
}

milo1.cs

Listing 16.19 milo1.cs Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    milo1.cs
//
//    Purpose:  Demonstrates how to solve a small mixed
//              integer linear optimization problem.
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class milo1
  {
    public static void Main(string[] args)
    {
      double[][] A =
      { new double[] { 50.0, 31.0 },
        new double[] { 3.0,  -2.0 }
      };
      double[] c = { 1.0, 0.64 };
      using (Model M = new Model("milo1"))
      {
        Variable x = M.Variable("x", 2, Domain.Integral(Domain.GreaterThan(0.0)));

        // Create the constraints
        //      50.0 x[0] + 31.0 x[1] <= 250.0
        //       3.0 x[0] -  2.0 x[1] >= -4.0
        M.Constraint("c1", Expr.Dot(A[0], x), Domain.LessThan(250.0));
        M.Constraint("c2", Expr.Dot(A[1], x), Domain.GreaterThan(-4.0));

        // Set max solution time
        M.SetSolverParam("mioMaxTime", 60.0);
        // Set max relative gap (to its default value)
        M.SetSolverParam("mioTolRelGap", 1e-4);
        // Set max absolute gap (to its default value)
        M.SetSolverParam("mioTolAbsGap", 0.0);

        // Set the objective function to (c^T * x)
        M.Objective("obj", ObjectiveSense.Maximize, Expr.Dot(c, x));

        // Solve the problem
        M.Solve();

        // Get the solution values
        double[] sol = x.Level();
        Console.WriteLine("x1,x2 = {0}, {1}", sol[0], sol[1]);
        double miorelgap = M.GetSolverDoubleInfo("mioObjRelGap");
        double mioabsgap = M.GetSolverDoubleInfo("mioObjAbsGap");
        Console.WriteLine("MIP rel gap = {0} ({0})", miorelgap, mioabsgap);
      }
    }
  }
}

mioinitsol.cs

Listing 16.20 mioinitsol.cs Click here to download.
//
//    Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//    File:    mioinitsol.cs
//
//    Purpose:  Demonstrates how to solve a small mixed
//             integer linear optimization problem
//             providing an initial feasible solution.
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class mioinitsol
  {
    public static void Main(string[] args)
    {
      double[] c = { 7.0, 10.0, 1.0, 5.0 };
      int n = 4;
      using (Model M = new Model("mioinitsol"))
      {

        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));
        x.Slice(0,3).MakeInteger();

        // Create the constraint
        M.Constraint(Expr.Sum(x), Domain.LessThan(2.5));

        // Set the objective function to (c^T * x)
        M.Objective("obj", ObjectiveSense.Maximize, Expr.Dot(c, x));

        // Assign values to integer variables.
        // We only set a slice of x     
        double[] init_sol = { 1, 1, 0 };
        x.Slice(0,3).SetLevel( init_sol );

        // Request constructing the solution from integer variable values
        M.SetSolverParam("mioConstructSol", "on");

        // Solve the problem
        M.Solve();

        // Get the solution values
        double[] sol = x.Level();
        Console.Write("x = [");
        for(int i=0;i<n;i++)
        {
          Console.Write("{0}, ",sol[i]);
        }
        Console.WriteLine("]");

        // Was the initial solution used?
        int constr = M.GetSolverIntInfo("mioConstructSolution");
        double constrVal = M.GetSolverDoubleInfo("mioConstructSolutionObj");
        Console.WriteLine("Construct solution utilization: " + constr);
        Console.WriteLine("Construct solution objective: " +  constrVal);
      }
    }
  }
}

nearestcorr.cs

Listing 16.21 nearestcorr.cs Click here to download.
/*
   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

   File:      nearestcorr.cs

   Purpose:
   Solves the nearest correlation matrix problem

     minimize   || A - X ||_F   s.t.  diag(X) = e, X is PSD

   as the equivalent conic program

     minimize     t

     subject to   (t, vec(A-X)) in Q
                  diag(X) = e
                  X >= 0
*/
using mosek.fusion;
using System;

namespace mosek.fusion.example
{
  public class nearestcorr
  {
    /** Assuming that e is an NxN expression, return the lower triangular part as a vector.
    */
    public static Expression Vec(Expression e)
    {
      int N       = e.GetShape()[0];
      int[] msubi = new int[N * (N + 1) / 2],
      msubj = new int[N * (N + 1) / 2];
      double[] mcof = new double[N * (N + 1) / 2];

      for (int i = 0, k = 0; i < N; ++i)
        for (int j = 0; j < i + 1; ++j, ++k)
        {
          msubi[k] = k;
          msubj[k] = i * N + j;
          if (i == j) mcof[k] = 1.0;
          else        mcof[k] = Math.Sqrt(2);
        }

      var S = Matrix.Sparse(N * (N + 1) / 2, N * N, msubi, msubj, mcof);
      return Expr.Mul(S, Expr.Flatten(e));
    }

    public static void Main(string[] argv)
    {
      int N = 5;
      var A = Matrix.Dense( new double[,]
      { { 0.0,  0.5,  - 0.1,  -0.2,   0.5},
        { 0.5,  1.25, -0.05, -0.1,   0.25},
        { -0.1, -0.05,  0.51,  0.02, -0.05},
        { -0.2, -0.1,   0.02,  0.54, -0.1},
        { 0.5,  0.25, -0.05, -0.1,   1.25}
      });

      var gammas = Range(0.0, 2.0, 0.1);
      var res  = new double[gammas.Length];
      var rank = new int[gammas.Length];

      nearestcorr_frobenius(A);
      nearestcorr_nn(A, gammas, res, rank);

      for (int i = 0; i < gammas.Length; ++i)
        Console.WriteLine("gamma = {0:f1}, res={1:e3}, rank = {2}", gammas[i], res[i], rank[i]);
    }

    public static void nearestcorr_nn(Matrix A, double[] gammas, double[] res, int[] rank)
    {
      int N = A.NumRows();
      using (var M = new Model("NucNorm"))
      {
        // Setup variables
        var t = M.Variable("t", 1, Domain.Unbounded());
        var X = M.Variable("X", Domain.InPSDCone(N));
        var w = M.Variable("w", N, Domain.GreaterThan(0.0));

        // (t, vec (X + diag(w) - A)) in Q
        var D = Expr.MulElm( Matrix.Eye(N), Var.Repeat(w, 1, N) );
        M.Constraint( Expr.Vstack( t, Vec(Expr.Sub(Expr.Add(X, D), A)) ), Domain.InQCone() );

        for (var k = 0; k < gammas.Length; ++k)
        {
          // Objective: Minimize t + gamma*Tr(X)
          var gamm_trX = Expr.Mul( gammas[k], Expr.Sum(X.Diag()) );
          M.Objective(ObjectiveSense.Minimize, Expr.Add(t, gamm_trX));
          M.Solve();

          // Find the eigenvalues of X and approximate rank
          var d = new double[N];
          mosek.LinAlg.syeig(mosek.uplo.lo, N, X.Level(), d);
          var rnk = 0; foreach (var v in d) if (v > 1e-6) ++rnk;

          res[k] = t.Level()[0];
          rank[k] = rnk;
        }
      }
    }

    public static void nearestcorr_frobenius(Matrix A)
    {
      int N = A.NumRows();
      using (var M = new Model("NearestCorrelation"))
      {
        // Setting up the variables
        var X = M.Variable("X", Domain.InPSDCone(N));
        var t = M.Variable("t", 1, Domain.Unbounded());

        // (t, vec (A-X)) \in Q
        M.Constraint( Expr.Vstack(t, Vec(Expr.Sub(A, X))), Domain.InQCone() );

        // diag(X) = e
        M.Constraint(X.Diag(), Domain.EqualsTo(1.0));

        // Objective: Minimize t
        M.Objective(ObjectiveSense.Minimize, t);

        // Solve the problem
        M.Solve();

        // Get the solution values
        Console.WriteLine("X = \n{0}", mattostr(X.Level(), N));
        Console.WriteLine("t = {0}", mattostr(t.Level(), N));
      }
    }

    /* Utility methods */
    private static string mattostr(double[] a, int n)
    {
      var c = new System.Globalization.CultureInfo("en-US");
      var b = new System.Text.StringBuilder();
      if (a.Length == 1 )
      {
        b.Append(String.Format(c, "{0:e3}", a[0]));
      }
      else
      {
        for (int i = 0; i < a.Length; ++i)
        {
          if (  i  % n == 0 )
            b.Append("[");
          b.Append(" ").Append(String.Format(c, "{0:e3}", a[i]));
          if ( ( i + 1 ) % n == 0 )
            b.Append("]\n");
        }
      }
      return b.ToString();
    }

    public static double[] Range (double start, double stop, double step)
    {
      int len;
      if (start < stop && step > 0.0)
        len = 1 + (int)((stop - start - 1) / step);
      else if (stop < start && step < 0)
        len = 1 + (int)((start - stop - 1) / (- step));
      else
        len = 0;
      double[] res = new double[len];
      double v = start;
      for (int i = 0; i < len; ++i, v += step)
        res[i] = v;
      return res;
    }
  }
}

opt_server_sync.cs

Listing 16.22 opt_server_sync.cs Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      opt_server_sync.cs
//
//  Purpose :   Demonstrates how to use MOSEK OptServer
//              to solve optimization problem synchronously
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class opt_server_sync
  {
    public static void Main(string[] args)
    {
  
      if (args.Length<1) {
        Console.WriteLine("Missing argument, syntax is:");
        Console.WriteLine("   opt_server_sync host [certpath]");
        return;
      }

      String serveraddr = args[0];
      String tlscert = (args.Length==2) ? args[1] : "";

      // Setup a simple test problem
      Model M = new Model("testOptServer");
      Variable x = M.Variable("x", 3, Domain.GreaterThan(0.0));
      M.Constraint("lc", Expr.Dot(new double[] {1.0, 1.0, 2.0}, x), Domain.EqualsTo(1.0));
      M.Objective("obj", ObjectiveSense.Minimize, Expr.Sum(x));

      // Attach log handler
      M.SetLogHandler(Console.Out);

      // Set OptServer URL
      M.OptserverHost(serveraddr);

      // Path to certificate, if any
      M.SetSolverParam("remoteTlsCertPath", tlscert);

      // Solve the problem on the OptServer
      M.Solve();

      // Get the solution
      double[] solx = x.Level();
      Console.WriteLine("x1,x2,x3 = {0}, {1}, {2}", solx[0], solx[1], solx[2]);
    }
  }
}

parallel.cs

Listing 16.23 parallel.cs Click here to download.
/*
  File : parallel.cs

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Purpose: Demonstrates parallel optimization using solveBatch()
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class Parallel
  {
    public static Model makeToyParameterizedModel() 
    {
      Model M = new Model();
      Variable x = M.Variable("x",3);
      Parameter p = M.Parameter("p");
      M.Objective(ObjectiveSense.Maximize, Expr.Sum(x));
      M.Constraint(Expr.Vstack(p,x), Domain.InQCone());
      return M;
    }


   /** Example of how to use Model.solveBatch()
    */
    public static void Main(string[] argv)
    {
      // Choose some sample parameters
      int n = 10;                 // Number of models to optimize
      int threadpoolsize = 4;     // Total number of threads available
      int threadspermodel = 1;    // Number of threads per each model

      // Create a toy model for this example
      Model M = makeToyParameterizedModel();

      // Set up n copies of the model with different data
      Model[] models = new Model[n];
      for(int i = 0; i < n ; i++)
      {
        models[i] = M.Clone();
        models[i].GetParameter("p").SetValue(i+1);
        // We can set the number of threads individually per model
        models[i].SetSolverParam("numThreads", threadspermodel);
      }

      // Solve all models in parallel
      SolverStatus[] status = Model.SolveBatch(false,         // No race
                                               -1.0,          // No time limit
                                               threadpoolsize,
                                               models);       // Array of Models to solve
      
      // Access the soutions
      for(int i = 0; i < n; i++) 
        if (status[i] == SolverStatus.OK)
          Console.WriteLine("Model  {0}: Status {1}  Solution Status {2}   objective  {3}  time {4}", 
            i, 
            status[i],
            models[i].GetPrimalSolutionStatus(),
            models[i].PrimalObjValue(),
            models[i].GetSolverDoubleInfo("optimizerTime"));
        else
          Console.WriteLine("Model  {0}: not solved", i);
    }
  }
}

parameters.cs

Listing 16.24 parameters.cs Click here to download.
/*
  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File :      parameters.cs

  Purpose :   Demonstrates a very simple example about how to set
              parameters and read information items
              with MOSEK Fusion
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class parameters
  {
    public static void Main(string[] args)
    {
      // Create the Model        
      using (Model M = new Model()) {
        Console.WriteLine("Test MOSEK parameter get/set functions");

        // Set log level (integer parameter)
        M.SetSolverParam("log", 1);
        // Select interior-point optimizer... (parameter with symbolic string values)
        M.SetSolverParam("optimizer", "intpnt");
        // ... without basis identification (parameter with symbolic string values)
        M.SetSolverParam("intpntBasis", "never");
        // Set relative gap tolerance (double parameter)
        M.SetSolverParam("intpntCoTolRelGap", 1.0e-7);

        // The same in a different way
        // M.SetSolverParam("intpntCoTolRelGap", "1.0e-7");

        // Incorrect value
        try {
          M.SetSolverParam("intpntCoTolRelGap", -1);
        }
        catch (mosek.fusion.ParameterError e) {
          Console.WriteLine("Wrong parameter value");
        }


        // Define and solve an optimization problem here
        // M.Solve()
        // After optimization: 

        Console.WriteLine("Get MOSEK information items");

        double tm = M.GetSolverDoubleInfo("optimizerTime");
        int it = M.GetSolverIntInfo("intpntIter");

        Console.WriteLine("Time: " + tm);
        Console.WriteLine("Iterations: " + it);
      }
    }
  }
}

portfolio_1_basic.cs

Listing 16.25 portfolio_1_basic.cs Click here to download.
/*
  File : portfolio_1_basic.cs

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model.
  
  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;

namespace mosek.fusion.example
{
  public class portfolio_1_basic
  {
    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i] * y[i];
      return r;
    }

    /*
    Purpose:
        Computes the optimal portfolio for a given risk

    Input:
        n: Number of assets
        mu: An n dimmensional vector of expected returns
        GT: A matrix with n columns so (GT")*GT  = covariance matrix"
        x0: Initial holdings
        w: Initial cash holding
        gamma: Maximum risk (=std. dev) accepted

    Output:
        Optimal expected return and the optimal portfolio
    */
    public static double BasicMarkowitz
    ( int n,
      double[] mu,
      double[,]GT,
      double[] x0,
      double   w,
      double   gamma)
    {

      using( Model M = new Model("Basic Markowitz"))
      {

        // Redirect log output from the solver to stdout for debugging.
        // if uncommented.
        //M.SetLogHandler(Console.Out);

        // Defines the variables (holdings). Shortselling is not allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));

        //  Maximize expected return
        M.Objective("obj", ObjectiveSense.Maximize, Expr.Dot(mu, x));

        // The amount invested  must be identical to intial wealth
        M.Constraint("budget", Expr.Sum(x), Domain.EqualsTo(w + sum(x0)));

        // Imposes a bound on the risk
        M.Constraint("risk", Expr.Vstack(gamma, Expr.Mul(GT, x)), Domain.InQCone());

        // Solves the model.
        M.Solve();

        return dot(mu, x.Level());
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {
      int        n      = 8;
      double     w      = 59.0;
      double[]   mu     = {0.07197349, 0.15518171, 0.17535435, 0.0898094 , 0.42895777, 0.39291844, 0.32170722, 0.18378628};
      double[]   x0     = {8.0, 5.0, 3.0, 5.0, 2.0, 9.0, 3.0, 6.0};
      double[]   gammas = {36};
      double[,]  GT     = {
        {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
        {0.0    , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
        {0.0    , 0.0    , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
        {0.0    , 0.0    , 0.0    , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.36096, 0.12574, 0.10157, 0.0571 },
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.21552, 0.05663, 0.06187},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.22514, 0.03327},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.2202 }
      };

      Console.WriteLine("\n-------------------------------------------------------------------");
      Console.WriteLine("Basic Markowitz portfolio optimization");
      Console.WriteLine("---------------------------------------------------------------------");
      foreach (var gamma in gammas)
      {
        double res = BasicMarkowitz(n, mu, GT, x0, w, gamma);
        Console.WriteLine("Expected return: {0,-12:f4}  St deviation: {1,-12:f4} ", res, gamma);
      }
    }
  }
}

portfolio_2_frontier.cs

Listing 16.26 portfolio_2_frontier.cs Click here to download.
/*
  File : portfolio_2_frontier.cs

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model.
                 Computes points on the efficient frontier.

  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;

namespace mosek.fusion.example
{
  public class portfolio_2_frontier
  {
    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i] * y[i];
      return r;
    }

    /*
      Purpose:
          Computes several portfolios on the optimal portfolios by

              for alpha in alphas:
                  maximize   expected return - alpha * variance
                  subject to the constraints

      Input:
          n: Number of assets
          mu: An n dimmensional vector of expected returns
          GT: A matrix with n columns so (GT")*GT  = covariance matrix"
          x0: Initial holdings
          w: Initial cash holding
          alphas: List of the alphas

      Output:
          The efficient frontier as list of tuples (alpha, expected return, variance)
     */
    public static void EfficientFrontier
    ( int n,
      double[]    mu,
      double [,]  GT,
      double[]    x0,
      double      w,
      double[]    alphas,
      double[]    frontier_mux,
      double[]    frontier_s)
    {
      using(Model M = new Model("Efficient frontier"))
      {
        //M.SetLogHandler(Console.Out);

        // Defines the variables (holdings). Shortselling is not allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0)); // Portfolio variables
        Variable s = M.Variable("s", 1, Domain.Unbounded());      // Variance variable

        M.Constraint("budget", Expr.Sum(x), Domain.EqualsTo(w + sum(x0)));

        // Computes the risk
        M.Constraint("variance", Expr.Vstack(s, 0.5, Expr.Mul(GT, x)), Domain.InRotatedQCone());

        //  Define objective as a weighted combination of return and variance
        Parameter alpha = M.Parameter();
        M.Objective("obj", ObjectiveSense.Maximize, Expr.Sub(Expr.Dot(mu, x), Expr.Mul(alpha, s)));

        // Solve the same problem for many values of parameter alpha
        for (int i = 0; i < alphas.Length; ++i) {
          alpha.SetValue(alphas[i]);
          M.Solve();

          frontier_mux[i] = dot(mu, x.Level());
          frontier_s[i]   = s.Level()[0];
        }
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {
      int        n      = 8;
      double     w      = 1.0;
      double[]   mu     = {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379};
      double[]   x0     = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
      double[,]  GT     = {
        {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
        {0.0    , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
        {0.0    , 0.0    , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
        {0.0    , 0.0    , 0.0    , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.36096, 0.12574, 0.10157, 0.0571 },
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.21552, 0.05663, 0.06187},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.22514, 0.03327},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.2202 }
      };

      // Some predefined alphas are chosen
      double[] alphas = { 0.0, 0.01, 0.1, 0.25, 0.30, 0.35, 0.4, 0.45, 0.5, 0.75, 1.0, 1.5, 2.0, 3.0, 10.0 };
      int      niter = alphas.Length;
      double[] frontier_mux = new double[niter];
      double[] frontier_s   = new double[niter];

      EfficientFrontier(n, mu, GT, x0, w, alphas, frontier_mux, frontier_s);
      Console.WriteLine("\n-------------------------------------------------------------------------");
      Console.WriteLine("Efficient frontier") ;
      Console.WriteLine("------------------------------------------------------------------------");
      Console.WriteLine("{0,-12}  {1,-12}  {2,-12}", "alpha", "return", "std. dev.") ;
      for (int i = 0; i < frontier_mux.Length; ++i)
        Console.WriteLine("{0,-12:f4}  {1,-12:e4}  {2,-12:e4}", alphas[i], frontier_mux[i], Math.Sqrt(frontier_s[i]));

    }
  }
}

portfolio_3_impact.cs

Listing 16.27 portfolio_3_impact.cs Click here to download.
/*
  File : portfolio_3_impact.cs

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model
                 with transaction costs of type x^(3/2)

  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;

namespace mosek.fusion.example
{
  public class portfolio_3_impact
  {
    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i] * y[i];
      return r;
    }


    /*
        Description:
            Extends the basic Markowitz model with a market cost term.

        Input:
            n: Number of assets
            mu: An n dimmensional vector of expected returns
            GT: A matrix with n columns so (GT')*GT  = covariance matrix'
            x0: Initial holdings
            w: Initial cash holding
            gamma: Maximum risk (=std. dev) accepted
            m: It is assumed that  market impact cost for the j'th asset is
               m_j|x_j-x0_j|^3/2

        Output:
           Optimal expected return and the optimal portfolio

    */
    public static void MarkowitzWithMarketImpact
    ( int n,
      double[] mu,
      double[,]GT,
      double[] x0,
      double   w,
      double   gamma,
      double[] m,
      double[] xsol,
      double[] tsol)
    {
      using(Model M = new Model("Markowitz portfolio with market impact"))
      {

        //M.SetLogHandler(Console.Out);

        // Defines the variables. No shortselling is allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));

        // Variables computing market impact
        Variable t = M.Variable("t", n, Domain.Unbounded());

        //  Maximize expected return
        M.Objective("obj", ObjectiveSense.Maximize, Expr.Dot(mu, x));

        // Invested amount + slippage cost = initial wealth
        M.Constraint("budget", Expr.Add(Expr.Sum(x), Expr.Dot(m, t)), Domain.EqualsTo(w + sum(x0)));

        // Imposes a bound on the risk
        M.Constraint("risk", Expr.Vstack(gamma, Expr.Mul(GT, x)), Domain.InQCone());

        // t >= |x-x0|^1.5 using a power cone
        M.Constraint("tz", Expr.Hstack(t, Expr.ConstTerm(n, 1.0), Expr.Sub(x,x0)), Domain.InPPowerCone(2.0/3.0));

        M.Solve();

        if (xsol != null)
          Array.Copy(x.Level(), xsol, n);
        if (tsol != null)
          Array.Copy(t.Level(), tsol, n);
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {
      int        n      = 8;
      double     w      = 1.0;
      double[]   mu     = {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379};
      double[]   x0     = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
      double[,]  GT     = {
        {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
        {0.0    , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
        {0.0    , 0.0    , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
        {0.0    , 0.0    , 0.0    , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.36096, 0.12574, 0.10157, 0.0571 },
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.21552, 0.05663, 0.06187},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.22514, 0.03327},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.2202 }
      };

      // Somewhat arbirtrary choice of m
      double[] m = new double[n]; for (int i = 0; i < n; ++i) m[i] = 0.01;
      double[] x = new double[n];
      double[] t = new double[n];
      double gamma = 0.36;

      MarkowitzWithMarketImpact(n, mu, GT, x0, w, gamma, m, x, t);
      Console.WriteLine("\n-----------------------------------------------------------------------");
      Console.WriteLine("Markowitz portfolio optimization with market impact cost");
      Console.WriteLine("------------------------------------------------------------------------");
      Console.WriteLine("Expected return: {0:e4} St deviation: {1:e4} Market impact cost: {2:e4}\n",
                        dot(mu, x),
                        gamma,
                        dot(m, t));

    }
  }
}

portfolio_4_transcost.cs

Listing 16.28 portfolio_4_transcost.cs Click here to download.
/*
  File : portfolio_4_transcost.cs

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model
                 with fixed setup costs and transaction costs
                 as a mixed-integer problem.

  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;

namespace mosek.fusion.example
{
  public class portfolio_4_transcost
  {
    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i] * y[i];
      return r;
    }


    /*
        Description:
            Extends the basic Markowitz model with a market cost term.

        Input:
            n: Number of assets
            mu: An n dimmensional vector of expected returns
            GT: A matrix with n columns so (GT")*GT  = covariance matrix"
            x0: Initial holdings
            w: Initial cash holding
            gamma: Maximum risk (=std. dev) accepted
            f: If asset j is traded then a fixed cost f_j must be paid
            g: If asset j is traded then a cost g_j must be paid for each unit traded

        Output:
           Optimal expected return and the optimal portfolio

    */
    public static double[] MarkowitzWithTransactionsCost
    ( int n,
      double[] mu,
      double[,]GT,
      double[] x0,
      double   w,
      double   gamma,
      double[] f,
      double[] g)
    {

      // Upper bound on the traded amount
      double[] u = new double[n];
      {
        double v = w + sum(x0);
        for (int i = 0; i < n; ++i) u[i] = v;
      }

      using( Model M = new Model("Markowitz portfolio with transaction costs") )
      {
        //M.SetLogHandler(Console.Out);

        // Defines the variables. No shortselling is allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));

        // Addtional "helper" variables
        Variable z = M.Variable("z", n, Domain.Unbounded());
        // Binary varables
        Variable y = M.Variable("y", n, Domain.Binary());

        //  Maximize expected return
        M.Objective("obj", ObjectiveSense.Maximize, Expr.Dot(mu, x));

        // Invest amount + transactions costs = initial wealth
        M.Constraint("budget", Expr.Add(Expr.Add(Expr.Sum(x), Expr.Dot(f, y)), Expr.Dot(g, z)),
                     Domain.EqualsTo(w + sum(x0)));

        // Imposes a bound on the risk
        M.Constraint("risk", Expr.Vstack(gamma, Expr.Mul(GT, x)), Domain.InQCone());

        // z >= |x-x0|
        M.Constraint("buy",  Expr.Sub(z, Expr.Sub(x, x0)), Domain.GreaterThan(0.0));
        M.Constraint("sell", Expr.Sub(z, Expr.Sub(x0, x)), Domain.GreaterThan(0.0));

        //M.constraint("trade", Expr.hstack(z,Expr.sub(x,x0)), Domain.inQcone())"

        // Consraints for turning y off and on. z-diag(u)*y<=0 i.e. z_j <= u_j*y_j
        M.Constraint("y_on_off", Expr.Sub(z, Expr.Mul(Matrix.Diag(u), y)), Domain.LessThan(0.0));

        // Integer optimization problems can be very hard to solve so limiting the
        // maximum amount of time is a valuable safe guard
        M.SetSolverParam("mioMaxTime", 180.0);
        M.Solve();
        return x.Level();
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {

      int        n      = 8;
      double     w      = 1.0;
      double[]   mu     = {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379};
      double[]   x0     = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
      double[,]  GT     = {
        {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
        {0.0    , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
        {0.0    , 0.0    , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
        {0.0    , 0.0    , 0.0    , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.36096, 0.12574, 0.10157, 0.0571 },
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.21552, 0.05663, 0.06187},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.22514, 0.03327},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.2202 }
    };

      double[] f = new double[n]; for (var i = 0; i < n; ++i) f[i] = 0.01;
      double[] g = new double[n]; for (var i = 0; i < n; ++i) g[i] = 0.001;
      double gamma = 0.36;

      double[] x = MarkowitzWithTransactionsCost(n, mu, GT, x0, w, gamma, f, g);
      Console.WriteLine("\n-------------------------------------------------------------------------");
      Console.WriteLine("Markowitz portfolio optimization with transaction cost");
      Console.WriteLine("------------------------------------------------------------------------");
      Console.WriteLine("Expected return: {0:e4}", dot(mu, x));

    }
  }
}

portfolio_5_card.cs

Listing 16.29 portfolio_5_card.cs Click here to download.
/*
  File : portfolio_5_card.cs

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Description :  Implements a basic portfolio optimization model
                 with cardinality constraints on number of assets traded.

  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;

namespace mosek.fusion.example
{
  public class portfolio_5_card
  {
    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i] * y[i];
      return r;
    }


    /*
        Description:
            Extends the basic Markowitz model with cardinality constraints.

        Input:
            n: Number of assets
            mu: An n dimmensional vector of expected returns
            GT: A matrix with n columns so (GT")*GT  = covariance matrix"
            x0: Initial holdings
            w: Initial cash holding
            gamma: Maximum risk (=std. dev) accepted
            k: Maximal number of assets in which we allow to change position.

        Output:
           Optimal expected return and the optimal portfolio

    */
    public static double[,] MarkowitzWithCardinality
    ( int n,
      double[] mu,
      double[,]GT,
      double[] x0,
      double   w,
      double   gamma,
      int[]    kValues)
    {
      // Upper bound on the traded amount
      double[] u = new double[n];
      {
        double v = w + sum(x0);
        for (int i = 0; i < n; ++i) u[i] = v;
      }

      using( Model M = new Model("Markowitz portfolio with cardinality bounds") )
      {
        //M.SetLogHandler(Console.Out);

        // Defines the variables. No shortselling is allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));

        // Addtional "helper" variables
        Variable z = M.Variable("z", n, Domain.Unbounded());
        // Binary varables
        Variable y = M.Variable("y", n, Domain.Binary());

        //  Maximize expected return
        M.Objective("obj", ObjectiveSense.Maximize, Expr.Dot(mu, x));

        // The amount invested  must be identical to initial wealth
        M.Constraint("budget", Expr.Sum(x), Domain.EqualsTo(w + sum(x0)));

        // Imposes a bound on the risk
        M.Constraint("risk", Expr.Vstack(gamma, Expr.Mul(GT, x)), Domain.InQCone());

        // z >= |x-x0|
        M.Constraint("buy",  Expr.Sub(z, Expr.Sub(x, x0)), Domain.GreaterThan(0.0));
        M.Constraint("sell", Expr.Sub(z, Expr.Sub(x0, x)), Domain.GreaterThan(0.0));

        // Consraints for turning y off and on. z-diag(u)*y<=0 i.e. z_j <= u_j*y_j
        M.Constraint("y_on_off", Expr.Sub(z, Expr.Mul(Matrix.Diag(u), y)), Domain.LessThan(0.0));

        // At most k assets change position
        Parameter cardMax = M.Parameter();
        M.Constraint("cardinality", Expr.Sub(Expr.Sum(y), cardMax), Domain.LessThan(0));

        // Integer optimization problems can be very hard to solve so limiting the
        // maximum amount of time is a valuable safe guard
        M.SetSolverParam("mioMaxTime", 180.0);
        
        // Solve multiple instances by varying the parameter k 
        double[,] results = new double[kValues.Length,n];

        for(int i = 0; i < kValues.Length; i++) {
          cardMax.SetValue(kValues[i]);
          M.Solve();
          double[] sol = x.Level();
          for(int j = 0; j < n; j++) results[i,j] = sol[j];
        }   

        return results;
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {
      int        n      = 8;
      double     w      = 1.0;
      double[]   mu     = {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379};
      double[]   x0     = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
      double[,]  GT     = {
        {0.30758, 0.12146, 0.11341, 0.11327, 0.17625, 0.11973, 0.10435, 0.10638},
        {0.0    , 0.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
        {0.0    , 0.0    , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
        {0.0    , 0.0    , 0.0    , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.36096, 0.12574, 0.10157, 0.0571 },
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.21552, 0.05663, 0.06187},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.22514, 0.03327},
        {0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.0    , 0.2202 }
      };
      double     gamma  = 0.25;

      int[]      kValues = { 1, 2, 3, 4, 5, 6, 7, 8 };

      Console.WriteLine("\n-------------------------------------------------------------------------");
      Console.WriteLine("Markowitz portfolio optimization with cardinality constraints");
      Console.WriteLine("------------------------------------------------------------------------");

      double[,] results = MarkowitzWithCardinality(n, mu, GT, x0, w, gamma, kValues);

      for(int K=1; K<=n; K++) 
      {
        Console.Write("Bound:  {0:d}  Solution: ", K);
        for(int i=0; i<n; i++)
        {
          Console.Write("{0:f4} ", results[K-1,i]);
        }
        Console.WriteLine();
      }
    }
  }
}

portfolio_6_factor.cs

Listing 16.30 portfolio_6_factor.cs Click here to download.
/*
  File : portfolio_6_factor.cs

  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  Purpose : Implements a portfolio optimization model using factor model.
  
  Note:
    This example uses LINQ, which is only available in .NET Framework 3.5 and later.
*/

using System.IO;
using System;
using System.Linq;
using System.Globalization;
using mosek;

namespace mosek.fusion.example
{
  public class portfolio_6_factor
  {
    public static double sum(double[] x)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i];
      return r;
    }

    public static double dot(double[] x, double[] y)
    {
      double r = 0.0;
      for (int i = 0; i < x.Length; ++i) r += x[i] * y[i];
      return r;
    }

    // Create diagonal matrix from an array.
    public static double[,] diag_matrix(double[] x) 
    {
      int n = x.Length;
      double[,] m = new double[n, n];  
      Array.Clear(m, 0, m.Length);
      for (int i = 0; i < n; ++i) m[i, i] = x[i];
      return m;
    }

    public static double[,] transpose(double[,] m) 
    {
      int ni = m.GetLength(0); 
      int nj = m.GetLength(1);
      double[,] mt = new double[nj, ni];  
      
      for (int i = 0; i < ni; ++i) 
      {
        for (int j = 0; j < nj; ++j) 
        {
          mt[j, i] = m[i, j];
        }
      }
      return mt;
    }

    public static double[,] array_sqrt(double[,] m) 
    {
      int ni = m.GetLength(0); 
      int nj = m.GetLength(1);
      double[,] sqrtm = new double[ni, nj];  
      
      for (int i = 0; i < ni; ++i) 
      {
        for (int j = 0; j < nj; ++j) 
        {
          sqrtm[i, j] = Math.Sqrt(m[i, j]);
        }
      }
      return sqrtm;
    }

    // Create horizontal block matrix
    public static double[,] concat_h(double[,] ml, double[,] mr) 
    {
      int ni = ml.GetLength(0);
      int njl = ml.GetLength(1);
      int njr = mr.GetLength(1);
      double[,] b = new double[ni, njl + njr];  
      
      for (int i = 0; i < ni; ++i) 
      {
        for (int j = 0; j < njl; ++j) b[i, j] = ml[i, j];
        for (int j = 0; j < njr; ++j) b[i, njl + j] = mr[i, j];
      }
      return b;
    }

    // Vectorize matrix (column-major order)
    public static double[] mat_to_vec_c(double[,] m) 
    {
      int ni = m.GetLength(0);
      int nj = m.GetLength(1);
      double[] c = new double[nj * ni];  
      
      for (int j = 0; j < nj; ++j) 
      {
        for (int i = 0; i < ni; ++i) 
        {
          c[j * ni + i] = m[i, j];
        }
      }
      return c;
    }

    // Reshape vector to matrix (column-major order)
    public static double[,] vec_to_mat_c(double[] c, int ni, int nj) 
    {
      double[,] m = new double[ni, nj];
      
      for (int j = 0; j < nj; ++j) 
      {
        for (int i = 0; i < ni; ++i) 
        {
          m[i, j] = c[j * ni + i];
        }
      }
      return m;
    }

    public static double[,] cholesky(double[,] m) 
    {
      int n = m.GetLength(0);
      double[] vecs = mat_to_vec_c(m);
      LinAlg.potrf(mosek.uplo.lo, n, vecs);
      double[,] s = vec_to_mat_c(vecs, n, n);

      // Zero out upper triangular part (LinAlg.Potrf does not use it, original matrix values remain there)
      for (int i = 0; i < n; ++i) 
      {
        for (int j = i+1; j < n; ++j) 
        {
          s[i, j] = 0.0;
        }
      }
      return s;
    }

    public static double[,] matrix_mul(double[,] a, double[,] b) 
    {
      int na = a.GetLength(0);
      int nb = b.GetLength(1);
      int k = b.GetLength(0);

      double[] vecm = new double[na * nb];
      Array.Clear(vecm, 0, vecm.Length);
      LinAlg.gemm(mosek.transpose.no, mosek.transpose.no, na, nb, k, 1.0, mat_to_vec_c(a), mat_to_vec_c(b), 1.0, vecm);
      double[,] m = vec_to_mat_c(vecm, na, nb);     
      
      return m;
    }

    /*
    Purpose:
        Computes the optimal portfolio for a given risk

    Input:
        n: Number of assets
        mu: An n dimmensional vector of expected returns
        GT: A matrix with n columns so (GT")*GT  = covariance matrix"
        x0: Initial holdings
        w: Initial cash holding
        gamma: Maximum risk (=std. dev) accepted

    Output:
        Optimal expected return and the optimal portfolio
    */
    public static double FactorMarkowitz
    ( int n,
      double[] mu,
      double[,]GT,
      double[] x0,
      double   w,
      double   gamma)
    {

      using( Model M = new Model("Factor Markowitz"))
      {

        // Redirect log output from the solver to stdout for debugging.
        // if uncommented.
        //M.SetLogHandler(Console.Out);

        // Defines the variables (holdings). Shortselling is not allowed.
        Variable x = M.Variable("x", n, Domain.GreaterThan(0.0));

        //  Maximize expected return
        M.Objective("obj", ObjectiveSense.Maximize, Expr.Dot(mu, x));

        // The amount invested  must be identical to intial wealth
        M.Constraint("budget", Expr.Sum(x), Domain.EqualsTo(w + sum(x0)));

        // Imposes a bound on the risk
        M.Constraint("risk", Expr.Vstack(gamma, Expr.Mul(GT, x)), Domain.InQCone());

        // Solves the model.
        M.Solve();

        return dot(mu, x.Level());
      }
    }


    /*
      The example. Reads in data and solves the portfolio models.
     */
    public static void Main(string[] argv)
    {
      int        n      = 8;
      double     w      = 1.0;
      double[]   mu     = {0.07197, 0.15518, 0.17535, 0.08981, 0.42896, 0.39292, 0.32171, 0.18379};
      double[]   x0     = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
      // Factor exposure matrix
      double[,]  B = 
      {
        {0.4256, 0.1869},
        {0.2413, 0.3877},
        {0.2235, 0.3697},
        {0.1503, 0.4612},
        {1.5325, -0.2633},
        {1.2741, -0.2613},
        {0.6939, 0.2372},
        {0.5425, 0.2116}
      };

      // Factor covariance matrix
      double[,]  S_F = 
      {
        {0.0620, 0.0577},
        {0.0577, 0.0908}
      };

      // Specific risk components
      double[]  resid_risk = {0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459};
      double[,] S_theta = diag_matrix(resid_risk);
      double[,] P = cholesky(S_F);
      double[,] BP = matrix_mul(B, P);  
      double[,] G = concat_h(BP, array_sqrt(S_theta));
      double[,] GT = transpose(G);
      double[]  gammas = {0.24, 0.28, 0.32, 0.36, 0.4, 0.44, 0.48};

      Console.WriteLine("\n-------------------------------------------------------------------");
      Console.WriteLine("Markowitz portfolio optimization with factor model");
      Console.WriteLine("---------------------------------------------------------------------");
      foreach (var gamma in gammas)
      {
        double res = FactorMarkowitz(n, mu, GT, x0, w, gamma);
        Console.WriteLine("Expected return: {0,-12:f4}  St deviation: {1,-12:f4} ", res, gamma);
      }
    }
  }
}

pow1.cs

Listing 16.31 pow1.cs Click here to download.
//
//   Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//   File:      pow1.cs
//
//   Purpose: Demonstrates how to solve the problem
//
//     maximize x^0.2*y^0.8 + z^0.4 - x
//           st x + y + 0.5z = 2
//              x,y,z >= 0
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class pow1
  {
    public static void Main(string[] args)
    {
      using (Model M = new Model("pow1"))
      {

        Variable x  = M.Variable("x", 3, Domain.Unbounded());
        Variable x3 = M.Variable();
        Variable x4 = M.Variable();

        // Create the linear constraint
        double[] aval = new double[] {1.0, 1.0, 0.5};
        M.Constraint(Expr.Dot(x, aval), Domain.EqualsTo(2.0));

        // Create the exponential conic constraint
        // Create the conic constraints
        M.Constraint(Var.Vstack(x.Slice(0,2), x3), Domain.InPPowerCone(0.2));
        M.Constraint(Expr.Vstack(x.Index(2), 1.0, x4), Domain.InPPowerCone(0.4));      

        // Set the objective function
        double[] cval = new double[] {1.0, 1.0, -1.0};
        M.Objective(ObjectiveSense.Maximize, Expr.Dot(cval, Var.Vstack(x3, x4, x.Index(0))));

        // Solve the problem
        M.Solve();

        // Get the linear solution values
        double[] solx = x.Level();
        Console.WriteLine("x,y,z = {0}, {1}, {2}", solx[0], solx[1], solx[2]);
      }
    }
  }
}

primal_svm.cs

Listing 16.32 primal_svm.cs Click here to download.
//
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      primal_svm.cs
//
//  Purpose: Implements a simple soft-margin SVM
//          using the Fusion API.

using System;
using System.Collections.Generic;

using mosek.fusion;

namespace mosek.fusion.example
{
  public class primal_svm
  {
    public static void Main(string[] args)
    {
      int nc = 10;
      int m = 50;
      int n = 3;
      int seed = 1;
      double variance = 1.0;
      double mean = 1.0;

      Random gen = new Random(seed);

      int nump = gen.Next(m);

      double [] y = new double[m];

      for (int i = 0; i < nump; i++)  y[i] = 1.0;
      for (int i = nump; i < m; i++) y[i] = -1.0;

      double [,] X = new double[m, n];

      for (int i = 0; i < nump; i++)
        for (int j = 0; j < n; j++) {
          double xx = Math.Sqrt( -2.0 * Math.Log(gen.NextDouble()) );
          double yy = Math.Sin( 2.0 * Math.PI * gen.NextDouble() );
          X[i, j] = (xx * yy) * variance + mean;
        }
      Model M = new Model("Primal SVM");
      try {
        Console.WriteLine("Number of data    : {0}\n", m);
        Console.WriteLine("Number of features: {0}\n", n);

        Variable w =  M.Variable( n, Domain.Unbounded() );
        Variable t =  M.Variable( 1, Domain.Unbounded() );
        Variable b =  M.Variable( 1, Domain.Unbounded() );
        Variable xi = M.Variable( m, Domain.GreaterThan(0.0) );

        M.Constraint( Expr.Add(Expr.MulElm( y, Expr.Sub( Expr.Mul(X, w), Var.Repeat(b, m) ) ), xi ),
                      Domain.GreaterThan( 1.0 ) );
        M.Constraint( Expr.Vstack(1.0, t, w) , Domain.InRotatedQCone());

        Console.WriteLine(" c     | b      | w");
        for (int i = 0; i < nc; i++) {
          double c = 500.0 * i;
          M.Objective(ObjectiveSense.Minimize, Expr.Add( t, Expr.Mul(c, Expr.Sum(xi) ) ) );
          M.Solve();

          Console.Write("{0} | {1} |", c, b.Level()[0] );
          for (int j = 0; j < n; j++)
            Console.Write(" {0}", w.Level()[j] );
          Console.WriteLine();
        }
      }
      finally {
        M.Dispose();
      }
    }
  }
}

qcqp_sdo_relaxation.cs

Listing 16.33 qcqp_sdo_relaxation.cs Click here to download.
/*
* Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
*
* File:      qcqp_sdo_relaxation.cs
*
* Purpose:   Demonstrate how to use SDP to solve
*            convex relaxation of a mixed-integer QCQP
**/
using System;

using mosek.fusion;
using mosek;

namespace mosek.fusion.example
{
  public class qcqp_sdo_relaxation
  {
    public static void Main(string[] args)
    {
      Random rnd = new Random();
      // problem dimensions
      int n = 20;
      int m = 2 * n;

      // problem data
      double[] A = new double[m * n];
      double[] b = new double[m];
      double[] c = new double[n];
      double[] P = new double[n * n];
      double[] q = new double[n];

      for (int j = 0; j < n; j++) {
        for (int i = 0; i < m; i++) {
          double u1 = rnd.NextDouble(); //these are uniform(0,1) random doubles
          double u2 = rnd.NextDouble();
          A[i * n + j] = Math.Sqrt(-2.0 * Math.Log(u1)) * Math.Sin(2.0 * Math.PI * u2);
        }
        c[j] = rnd.NextDouble();
      }

      // P = A^T A
      LinAlg.syrk(mosek.uplo.lo, mosek.transpose.yes,
                  n, m, 1.0, A, 0.0, P);
      for (int j = 0; j < n; j++) for (int i = j + 1; i < n; i++) P[i * n + j] = P[j * n + i];

      // q = -P c, b = A c
      LinAlg.gemv(mosek.transpose.no, n, n, -1.0, P, c, 0.0, q);
      LinAlg.gemv(mosek.transpose.no, m, n,  1.0, A, c, 0.0, b);

      // Solve the problems
      Model M = miqcqp_sdo_relaxation(n, Matrix.Dense(n, n, P), q);
      Model Mint = int_least_squares(n, Matrix.Dense(n, m, A).Transpose(), b);
      M.Solve();
      Mint.Solve();

      // rounded and optimal solution
      double[] xRound = M.GetVariable("Z").Slice(new int[] {0, n}, new int[] {n, n + 1}).Level();
      for (int i = 0; i < n; i++) xRound[i] = Math.Round(xRound[i]);
      double[] yRound = (double[]) b.Clone();
      double[] xOpt   = Mint.GetVariable("x").Level();
      double[] yOpt   = (double[]) b.Clone();
      LinAlg.gemv(mosek.transpose.no, m, n,  1.0, A, xRound, -1.0, yRound);        // Ax_round-b
      LinAlg.gemv(mosek.transpose.no, m, n,  1.0, A, xOpt, -1.0, yOpt);            // Ax_opt-b

      Console.WriteLine(M.GetSolverDoubleInfo("optimizerTime") + " " + Mint.GetSolverDoubleInfo("optimizerTime"));
      Console.WriteLine(Math.Sqrt(LinAlg.dot(m, yRound, yRound)) + " " +
                        Math.Sqrt(LinAlg.dot(m, yOpt, yOpt)));
    }

    public static Model miqcqp_sdo_relaxation(int n, Matrix P, double[] q) {
      Model M = new Model();

      Variable Z = M.Variable("Z", Domain.InPSDCone(n + 1));

      Variable X = Z.Slice(new int[] {0, 0}, new int[] {n, n});
      Variable x = Z.Slice(new int[] {0, n}, new int[] {n, n + 1});

      M.Constraint( Expr.Sub(X.Diag(), x), Domain.GreaterThan(0.0) );
      M.Constraint( Z.Index(n, n), Domain.EqualsTo(1.0) );

      M.Objective( ObjectiveSense.Minimize, Expr.Add(
                     Expr.Sum( Expr.MulElm( P, X ) ),
                     Expr.Mul( 2.0, Expr.Dot(x, q) )
                   ) );

      return M;
    }

    // A direct integer model for minimizing |Ax-b|
    public static Model int_least_squares(int n, Matrix A, double[] b) {
      Model M = new Model();

      Variable x = M.Variable("x", n, Domain.Integral(Domain.Unbounded()));
      Variable t = M.Variable("t", 1, Domain.Unbounded());

      M.Constraint( Expr.Vstack(t, Expr.Sub(Expr.Mul(A, x), b)), Domain.InQCone() );
      M.Objective( ObjectiveSense.Minimize, t );

      return M;
    }
  }
}

reoptimization.cs

Listing 16.34 reoptimization.cs Click here to download.
/**
  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File :      reoptimization.cs

  Purpose:   Demonstrates how to solve a  linear
             optimization problem using the MOSEK API
             and modify and re-optimize the problem.
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class reoptimization
  {
    public static void printsol(double[] x) {
      Console.Write("x = ");
      foreach(double val in x)
          Console.Write(String.Format("{0:F4} ", val));
      Console.WriteLine();
    }

    public static void Main(string[] args)
    {
      using (Model M = new Model())
      {
        double[]   c = new double[] { 1.5, 2.5, 3.0 };
        double[,]  A = new double[,] { {2, 4, 3},
                                       {3, 2, 3},
                                       {2, 3, 2} };
        double[]   b = new double[] { 100000.0, 50000.0, 60000.0 };
        int   numvar = c.Length;
        int   numcon = b.Length;

        // Create a model and input data
        Variable x = M.Variable(numvar, Domain.GreaterThan(0.0));
        Constraint con = M.Constraint(Expr.Mul(A, x), Domain.LessThan(b));
        M.Objective(ObjectiveSense.Maximize, Expr.Dot(c, x));
        // Solve the problem
        M.Solve();
        printsol(x.Level());

        /************** Change an element of the A matrix ****************/
        con.Index(0).Update(Expr.Mul(3.0, x.Index(0)), x.Index(0));
        M.Solve();
        printsol(x.Level());

        /*************** Add a new variable ******************************/
        // Create a variable and a compound view of all variables
        Variable x3 = M.Variable(Domain.GreaterThan(0.0));
        Variable xNew = Var.Vstack(x, x3);
        // Add to the exising constraint
        con.Update(Expr.Mul(x3, new double[]{4, 0, 1}), x3);
        // Change the objective to include x3
        M.Objective(ObjectiveSense.Maximize, Expr.Dot(new double[]{1.5,2.5,3.0,1.0}, xNew));
        M.Solve();
        printsol(xNew.Level());

        /**************** Add a new constraint *****************************/
        Constraint con2 = M.Constraint(Expr.Dot(xNew, new double[]{1, 2, 1, 1}), Domain.LessThan(30000.0));
        M.Solve();
        printsol(xNew.Level());

        /**************** Change constraint bounds *****************************/
        // Assemble all constraints we previously defined into one
        Constraint cAll = Constraint.Vstack(con, con2);
        // Change bounds by effectively updating fixed terms with the difference
        cAll.Update(new double[]{20000, 10000, 10000, 8000});
        M.Solve();
        printsol(xNew.Level());
      }
    }
  }
}

response.cs

Listing 16.35 response.cs Click here to download.
/*
  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File:      response.cs

  Purpose:   This example demonstrates proper response handling
             for problems solved with the interior-point optimizers.
*/
using System;
using mosek;
using mosek.fusion;
using System.Text;

namespace mosek.fusion.example
{
  public class response
  {
    // Set up a small artificial conic problem to test with
    public static void SetupExample(Model M) {
      Variable x = M.Variable("x", 3, Domain.GreaterThan(0.0));
      Variable y = M.Variable("y", 3, Domain.Unbounded());
      Variable z1 = Var.Vstack(y.Index(0), x.Slice(0, 2));
      Variable z2 = Var.Vstack(y.Slice(1, 3), x.Index(2));
      M.Constraint("lc", Expr.Dot(new double[] {1.0, 1.0, 2.0}, x), Domain.EqualsTo(1.0));
      M.Constraint("qc1", z1, Domain.InQCone());
      M.Constraint("qc2", z2, Domain.InRotatedQCone());
      M.Objective("obj", ObjectiveSense.Minimize, Expr.Sum(y));
    }

    public static void Main(String[] argv) {
      Model M = new Model();

      // (Optional) set a log stream
      // M.SetLogHandler(Console.Out);
      
      // (Optional) uncomment to see what happens when solution status is unknown
      // M.SetSolverParam("intpntMaxIterations", 1);

      // In this example we set up a small conic problem
      SetupExample(M);

      // Optimize
      try 
      {
        M.Solve();

        // We expect solution status OPTIMAL (this is also default)
        M.AcceptedSolutionStatus(AccSolutionStatus.Optimal);

        Variable x = M.GetVariable("x");
        long xsize = x.GetSize();
        double[] xVal = x.Level();
        Console.Write("Optimal value of x = ");
        for(int i = 0; i < xsize; ++i)
          Console.Write(xVal[i] + " ");
        Console.WriteLine("\nOptimal primal objective: " + M.PrimalObjValue());
        // .. continue analyzing the solution

      }
      catch (OptimizeError e)
      {
        Console.WriteLine("Optimization failed. Error: " + e.ToString());
      }
      catch (SolutionError)
      {
        // The solution with at least the expected status was not available.
        // We try to diagnoze why.
        Console.WriteLine("Requested solution was not available.");
        ProblemStatus prosta = M.GetProblemStatus();
        switch(prosta)
        {
          case ProblemStatus.DualInfeasible:
            Console.WriteLine("Dual infeasibility certificate found.");
            break;

          case ProblemStatus.PrimalInfeasible:
            Console.WriteLine("Primal infeasibility certificate found.");
            break;

          case ProblemStatus.Unknown:
            // The solutions status is unknown. The termination code
            // indicates why the optimizer terminated prematurely.
            Console.WriteLine("The solution status is unknown.");
            StringBuilder symname = new StringBuilder();
            StringBuilder desc = new StringBuilder();
            Env.getcodedesc((rescode)M.GetSolverIntInfo("optimizeResponse"), symname, desc);
            Console.WriteLine("  Termination code: {0} {1}", symname, desc);
            break;

          default:
            Console.WriteLine("Another unexpected problem status: " + prosta);
            break;
        }
      }
      catch (Exception e)
      {
        Console.WriteLine("Unexpected error: " + e.ToString());
      }

      M.Dispose();
    }
  }
}

sdo1.cs

Listing 16.36 sdo1.cs Click here to download.
//
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      sdo1.cs
//
//  Purpose:
//  Solves the mixed semidefinite and conic quadratic optimization problem
//
//                   [2, 1, 0]
//    minimize    Tr [1, 2, 1] * X + x0
//                   [0, 1, 2]
//
//                   [1, 0, 0]
//    subject to  Tr [0, 1, 0] * X + x0           = 1.0
//                   [0, 0, 1]
//
//                   [1, 1, 1]
//                Tr [1, 1, 1] * X      + x1 + x2 = 0.5
//                   [1, 1, 1]
//
//                   X >> 0,  x0 >= (x1^2 + x2^2) ^ (1/2)
//
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class sdo1
  {
    public static void Main(string[] args)
    {
      using (Model M  = new Model("sdo1"))
      {
        // Setting up the variables
        Variable X  = M.Variable("X", Domain.InPSDCone(3));
        Variable x  = M.Variable("x", Domain.InQCone(3));

        // Setting up constant coefficient matrices
        Matrix C  = Matrix.Dense ( new double[,] { {2, 1, 0}, {1, 2, 1}, {0, 1, 2}} );
        Matrix A1 = Matrix.Eye(3);
        Matrix A2 = Matrix.Ones(3,3);

        // Objective
        M.Objective(ObjectiveSense.Minimize, Expr.Add(Expr.Dot(C, X), x.Index(0)));

        // Constraints
        M.Constraint("c1", Expr.Add(Expr.Dot(A1, X), x.Index(0)), Domain.EqualsTo(1.0));
        M.Constraint("c2", Expr.Add(Expr.Dot(A2, X), Expr.Sum(x.Slice(1, 3))),
                     Domain.EqualsTo(0.5));

        M.Solve();

        Console.WriteLine("[{0}]", (new Utils.StringBuffer()).A(X.Level()).ToString());
        Console.WriteLine("[{0}]", (new Utils.StringBuffer()).A(x.Level()).ToString());
      }
    }
  }
}

sdo2.cs

Listing 16.37 sdo2.cs Click here to download.
/*
   Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
 
   File :      sdo2.cs
 
   Purpose :   Solves the semidefinite problem with two symmetric variables:
 
                  min   <C1,X1> + <C2,X2>
                  st.   <A1,X1> + <A2,X2> = b
                              (X2)_{1,2} <= k
                 
                  where X1, X2 are symmetric positive semidefinite,
 
                  C1, C2, A1, A2 are assumed to be constant symmetric matrices,
                  and b, k are constants.
*/
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class sdo2
  {
    public static void Main(string[] args)
    {

      // Sample data in sparse, symmetric triplet format
      int[]    C1_k = {0, 2};
      int[]    C1_l = {0, 2};
      double[] C1_v = {1, 6};
      int[]    A1_k = {0, 2, 0, 2};
      int[]    A1_l = {0, 0, 2, 2};
      double[] A1_v = {1, 1, 1, 2};
      int[]    C2_k = {0, 1, 0, 1, 2};
      int[]    C2_l = {0, 0, 1, 1, 2};
      double[] C2_v = {1, -3, -3, 2, 1};
      int[]    A2_k = {1, 0, 1, 3};
      int[]    A2_l = {0, 1, 1, 3};
      double[] A2_v = {1, 1, -1, -3};
      double b = 23;
      double k = -3;

      // Convert input data into Fusion sparse matrices
      Matrix C1 = Matrix.Sparse(3, 3, C1_k, C1_l, C1_v);
      Matrix C2 = Matrix.Sparse(4, 4, C2_k, C2_l, C2_v);
      Matrix A1 = Matrix.Sparse(3, 3, A1_k, A1_l, A1_v);
      Matrix A2 = Matrix.Sparse(4, 4, A2_k, A2_l, A2_v);

      using (Model M  = new Model("sdo2"))
      {
        // Two semidefinite variables
        Variable X1 = M.Variable(Domain.InPSDCone(3));
        Variable X2 = M.Variable(Domain.InPSDCone(4));

        // Objective
        M.Objective(ObjectiveSense.Minimize, Expr.Add(Expr.Dot(C1,X1), Expr.Dot(C2,X2)));

        // Equality constraint
        M.Constraint(Expr.Add(Expr.Dot(A1,X1), Expr.Dot(A2,X2)), Domain.EqualsTo(b));

        // Inequality constraint
        M.Constraint(X2.Index(new int[] {0,1}), Domain.LessThan(k));

        // Solve
        M.SetLogHandler(Console.Out);
        M.Solve();

        // Print solution
        Console.WriteLine("Solution (vectorized):");
        Console.WriteLine("[{0}]", (new Utils.StringBuffer()).A(X1.Level()).ToString());
        Console.WriteLine("[{0}]", (new Utils.StringBuffer()).A(X2.Level()).ToString());
      }
    }
  }
}

sdo3.cs

Listing 16.38 sdo3.cs Click here to download.
/*
  Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.

  File :      sdo3.cs

  Purpose :   Solves the semidefinite problem:

                 min   tr(X_1) + ... + tr(X_n)
                 st.   <A_11,X_1> + ... + <A_1n,X_n> >= b_1
                       ...
                       <A_k1,X_1> + ... + <A_kn,X_n> >= b_k
                
                 where X_i are symmetric positive semidefinite of dimension d,

                 A_ji are constant symmetric matrices and b_i are constant.

              This example is to demonstrate creating and using 
              many matrix variables of the same dimension.
*/

using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class sdo3
  {

    // A helper method computing a semidefinite slice of a 3-dim variable
    public static Variable Slice(Variable X, int d, int j) {
      return
        X.Slice(new int[] {j,0,0}, new int[] {j+1,d,d})
         .Reshape(new int[] {d,d});
    }

    public static void Main(string[] args)
    {

      // Sample input data
      int n = 100;
      int d = 4;
      int k = 3;
      double[] b = {9,10,11};
      double[][,] A = new double[n*k][,];
      var rand = new Random();
      for(int i=0; i<n*k; i++) {
        A[i] = new double[d,d];
        for(int s1=0; s1<d; s1++)
          for(int s2=0; s2<=s1; s2++)
            A[i][s1,s2] = A[i][s2,s1] = rand.NextDouble();
      }

      // Create a model with n semidefinite variables od dimension d x d
      using (Model M  = new Model("sdo3"))
      {
        Variable X = M.Variable(Domain.InPSDCone(d, n));

        // Pick indexes of diagonal entries for the objective
        int[,] alldiag = new int[d*n,3];
        for(int j=0; j<n; j++) for(int s=0; s<d; s++) {
          alldiag[j*d+s,0] = j; 
          alldiag[j*d+s,1] = alldiag[j*d+s,2] = s;
        }
        M.Objective(ObjectiveSense.Minimize, Expr.Sum( X.Pick(alldiag) ));

        // Each constraint is a sum of inner products
        // Each semidefinite variable is a slice of X
        for(int i=0; i< k; i++) {
          Expression[] addlist = new Expression[n];
          for(int j=0; j<n; j++) 
            addlist[j] = Expr.Dot(A[i*n+j], Slice(X, d, j));
          M.Constraint(Expr.Add(addlist), Domain.GreaterThan(b[i]));
        }

        // Solve
        M.SetLogHandler(Console.Out);           // Add logging
        M.WriteTask("sdo3.ptf");                // Save problem in readable format
        M.Solve();

        // Get results. Each variable is a slice of X
        Console.WriteLine("Contributing variables:");
        for(int j=0; j<n; j++) {
          double[] Xj = Slice(X, d, j).Level();
          double maxval = 0;
          for(int s=0; s<d*d; s++) maxval = Math.Max(maxval, Xj[s]);
          if (maxval > 1e-6) {
            Console.WriteLine("X" + j + "=");
            for(int s1=0; s1<d; s1++) {
              for(int s2=0; s2<d; s2++)
                Console.Write(Xj[s1*d+s1] + "  ");    
              Console.WriteLine();
            }
          }
        }
      } 
    }
  }
}

sospoly.cs

Listing 16.39 sospoly.cs Click here to download.
////
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      sospoly.cs
//
//  Purpose:
//  Models the cone of nonnegative polynomials and nonnegative trigonometric
//  polynomials using Nesterov's framework  [1].
//
//  Given a set of coefficients (x0, x1, ..., xn), the functions model the
//  cone of nonnegative polynomials
//
//  P_m = { x \in R^{n+1} | x0 + x1*t + ... xn*t^n >= 0, forall t \in I }
//
//  where I can be the entire real axis, the semi-infinite interval [0,inf), or
//  a finite interval I = [a, b], respectively.
//
//  References:
//
//  [1] "Squared Functional Systems and Optimization Problems",
//      Y. Nesterov, in High Performance Optimization,
//      Kluwer Academic Publishers, 2000.
using mosek.fusion;
using System;

namespace mosek.fusion.example
{
  class sospoly
  {
    // Creates a Hankel matrix of dimension n+1, where
    // H_lk = a if l+k=i, and 0 otherwise.
    public static Matrix Hankel(int n, int i, double a)
    {
      if (i < n + 1)
        return Matrix.Sparse(n + 1, n + 1, Range(i, -1, -1), Range(i + 1), a);
      else
        return Matrix.Sparse(n + 1, n + 1, Range(n, i - n - 1, -1), Range(i - n, n + 1), a);
    }

    // Models the cone of nonnegative polynomials on the real axis
    public static void nn_inf(Model M, Variable x)
    {
      int m = (int)x.GetSize() - 1;
      int n = (m / 2); // degree of polynomial is 2*n

      //assert(m == 2*n)

      // Setup variables
      Variable X = M.Variable(Domain.InPSDCone(n + 1));

      // x_i = Tr H(n, i) * X  i=0,...,m
      for (int i = 0; i < m + 1; ++i)
        M.Constraint( Expr.Sub(x.Index(i), Expr.Dot(Hankel(n, i, 1.0), X)), Domain.EqualsTo(0.0));
    }

    // Models the cone of nonnegative polynomials on the semi-infinite interval [0,inf)
    public static void nn_semiinf(Model M, Variable x)
    {
      int n = (int)x.GetSize() - 1;
      int n1 = n / 2,
          n2 = (n - 1) / 2;

      // Setup variables
      Variable X1 = M.Variable(Domain.InPSDCone(n1 + 1));
      Variable X2 = M.Variable(Domain.InPSDCone(n2 + 1));

      // x_i = Tr H(n1, i) * X1 + Tr H(n2,i-1) * X2, i=0,...,n
      for (int i = 0; i < n + 1; ++i)
        M.Constraint( Expr.Sub(x.Index(i),
                               Expr.Add(Expr.Dot(Hankel(n1, i, 1.0),  X1),
                                        Expr.Dot(Hankel(n2, i - 1, 1.0), X2))), Domain.EqualsTo(0.0) );
    }

    // Models the cone of nonnegative polynomials on the finite interval [a,b]
    public static void nn_finite(Model M, Variable x, double a, double b)
    {
      //assert(a < b)
      int m = (int)x.GetSize() - 1;
      int n = m / 2;

      if (m == 2 * n)
      {
        Variable X1 = M.Variable(Domain.InPSDCone(n + 1));
        Variable X2 = M.Variable(Domain.InPSDCone(n));

        // x_i = Tr H(n,i)*X1 + (a+b)*Tr H(n-1,i-1) * X2 - a*b*Tr H(n-1,i)*X2 - Tr H(n-1,i-2)*X2, i=0,...,m
        for (int i = 0; i < m + 1; ++i)
          M.Constraint( Expr.Sub(x.Index(i),
                                 Expr.Add(Expr.Sub(Expr.Dot(Hankel(n, i, 1.0),      X1), Expr.Dot(Hankel(n - 1, i, a * b), X2)),
                                          Expr.Sub(Expr.Dot(Hankel(n - 1, i - 1, a + b), X2), Expr.Dot(Hankel(n - 1, i - 2, 1.0),  X2)))),
                        Domain.EqualsTo(0.0) );
      }
      else
      {
        Variable X1 = M.Variable(Domain.InPSDCone(n + 1));
        Variable X2 = M.Variable(Domain.InPSDCone(n + 1));

        // x_i = Tr H(n,i-1)*X1 - a*Tr H(n,i)*X1 + b*Tr H(n,i)*X2 - Tr H(n,i-1)*X2, i=0,...,m
        for (int i = 0; i < m + 1; ++i)
          M.Constraint( Expr.Sub(x.Index(i),
                                 Expr.Add(Expr.Sub(Expr.Dot(Hankel(n, i - 1, 1.0),  X1), Expr.Dot(Hankel(n, i, a), X1)),
                                          Expr.Sub(Expr.Dot(Hankel(n, i, b),     X2), Expr.Dot(Hankel(n, i - 1, 1.0),  X2)))),
                        Domain.EqualsTo(0.0) );
      }
    }

    // returns variables u representing the derivative of
    //  x(0) + x(1)*t + ... + x(n)*t^n,
    // with u(0) = x(1), u(1) = 2*x(2), ..., u(n-1) = n*x(n).
    public static Variable diff(Model M, Variable x)
    {
      int n = (int)x.GetSize() - 1;
      Variable u = M.Variable(n, Domain.Unbounded());
      M.Constraint(Expr.Sub(u, Expr.MulElm(Matrix.Dense(n, 1, Range(1.0, n + 1)), x.Slice(1, n + 1))), Domain.EqualsTo(0.0));
      return u;
    }

    public static double[] fitpoly(double[,] data, int n)
    {
      using (var M = new Model("smooth poly"))
      {
        int m = data.GetLength(0);
        double[] Adata = new double[m * (n + 1)];
        for (int j = 0, k = 0; j < m; ++j)
          for (int i = 0; i < n + 1; ++i, ++k)
            Adata[k] = Math.Pow(data[j, 0], i);

        Matrix A = Matrix.Dense(m, n + 1, Adata);
        double[] b = new double[m];
        for (int j = 0; j < m; ++j)
          b[j] = data[j, 1];

        Variable x  = M.Variable("x", n + 1, Domain.Unbounded());
        Variable z  = M.Variable("z", 1,   Domain.Unbounded());
        Variable dx = diff(M, x);

        M.Constraint(Expr.Mul(A, x), Domain.EqualsTo(b));

        // z - f'(t) >= 0, for all t \in [a, b]
        Variable ub = M.Variable(n, Domain.Unbounded());
        M.Constraint(Expr.Sub(ub,
                              Expr.Vstack(Expr.Sub(z, dx.Index(0)), Expr.Neg(dx.Slice(1, n)))),
                     Domain.EqualsTo(0.0));
        nn_finite(M, ub, data[0, 0], data[m - 1, 0]);

        // f'(t) + z >= 0, for all t \in [a, b]
        Variable lb = M.Variable(n, Domain.Unbounded());
        M.Constraint(Expr.Sub(lb,
                              Expr.Vstack(Expr.Add(z, dx.Index(0)), dx.Slice(1, n).AsExpr())),
                     Domain.EqualsTo(0.0));
        nn_finite(M, lb, data[0, 0], data[m - 1, 0]);

        M.Objective(ObjectiveSense.Minimize, z);
        M.Solve();
        return x.Level();
      }
    }
    public static void Main(string[] argv)
    {

      double[,] data =
      { { -1.0, 1.0 },
        {  0.0, 0.0 },
        {  1.0, 1.0 }
      };

      double[] x2 = fitpoly(data, 2);
      double[] x4 = fitpoly(data, 4);
      double[] x8 = fitpoly(data, 8);

      Console.WriteLine("fitpoly(data,2) -> {0}", arrtostr(x2));
      Console.WriteLine("fitpoly(data,4) -> {0}", arrtostr(x4));
      Console.WriteLine("fitpoly(data,8) -> {0}", arrtostr(x8));
    }



    // Some utility functions to make things look nicer
    public static int[] Range (int start, int stop, int step)
    {
      int len;
      if (start < stop && step > 0)
        len = 1 + (stop - start - 1) / step;
      else if (stop < start && step < 0)
        len = 1 + (start - stop - 1) / (- step);
      else
        len = 0;
      int[] res = new int[len];
      for (int i = 0, v = start; i < len; ++i, v += step)
        res[i] = v;
      return res;
    }
    public static int[] Range (int start, int stop) { return Range(start, stop, 1); }
    public static int[] Range (int stop) { return Range(0, stop, 1); }
    public static double[] Range (double start, double stop, double step)
    {
      int len;
      if (start < stop && step > 0.0)
        len = 1 + (int)((stop - start - 1) / step);
      else if (stop < start && step < 0)
        len = 1 + (int)((start - stop - 1) / (- step));
      else
        len = 0;
      double[] res = new double[len];
      double v = start;
      for (int i = 0; i < len; ++i, v += step)
        res[i] = v;
      return res;
    }
    public static double[] Range (double start, double stop) { return Range(start, stop, 1.0); }
    public static double[] Range (double stop) { return Range(0.0, stop, 1.0); }

    public static double[] Fill (int num, double v) { double[] r = new double[num]; for (int i = 0; i < num; ++i) r[i] = v; return r;}

    private static string arrtostr(double[] a)
    {
      var c = new System.Globalization.CultureInfo("en-US");
      var b = new System.Text.StringBuilder("[");
      if (a.Length > 0) b.Append(String.Format(c, "{0:e3}", a[0]));
      for (int i = 1; i < a.Length; ++i) b.Append(",").Append(String.Format(c, "{0:e3}", a[i]));
      b.Append("]");
      return b.ToString();
    }
  }
}

sudoku.cs

Listing 16.40 sudoku.cs Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      sudoku.cs
//
// Purpose:  A MILP-based SUDOKU solver
//
using System;
using System.Collections.Generic;

using mosek.fusion;

namespace mosek.fusion.example
{
  public class sudoku
  {
    public static void Main(string[] args)
    {
      int m = 3;
      int n = m * m;

      //fixed cells in human readable (i.e. 1-based) format
      int[,] hr_fixed = {
        {1, 5, 4},
        {2, 2, 5}, {2, 3, 8}, {2, 6, 3},
        {3, 2, 1}, {3, 4, 2}, {3, 5, 8}, {3, 7, 9},
        {4, 2, 7}, {4, 3, 3}, {4, 4, 1}, {4, 7, 8}, {4, 8, 4},
        {6, 2, 4}, {6, 3, 1}, {6, 6, 9}, {6, 7, 2}, {6, 8, 7},
        {7, 3, 4}, {7, 5, 6}, {7, 6, 5}, {7, 8, 8},
        {8, 4, 4}, {8, 7, 1}, {8, 8, 6},
        {9, 5, 9}
      };
      int nf = hr_fixed.Length / 3;
      int[,] fixed_cells = new int[nf, 3];
      for (int i = 0; i < nf; i++)
        for (int d = 0; d < m; d++)
          fixed_cells[i, d] = hr_fixed[i, d] - 1;

      using (Model M = new Model("SUDOKU")) {
        M.SetLogHandler(Console.Out);
        Variable x = M.Variable(new int[] { n, n, n }, Domain.Binary());

        //each value only once per dimension
        for (int d = 0; d < m; d++)
          M.Constraint(Expr.Sum(x, d), Domain.EqualsTo(1.0));

        //each number must appears only once in a block
        for (int k = 0; k < n; k++)
          for (int i = 0; i < m; i++)
            for (int j = 0; j < m; j++)
            {
              Variable block = x.Slice(new int[] { i * m, j * m, k },
                                       new int[] { (i + 1)*m, (j + 1)*m, k + 1 });
              M.Constraint(Expr.Sum(block), Domain.EqualsTo(1.0));
            }

        M.Constraint(x.Pick(fixed_cells), Domain.EqualsTo(1.0));

        M.Solve();

        //print the solution, if any...
        if (M.GetPrimalSolutionStatus() == SolutionStatus.Optimal)
          print_solution(m, x);
        else
          Console.WriteLine("No solution found!\n");
      }
    }

    public static void print_solution(int m, Variable x)
    {
      int n = m * m;
      Console.WriteLine("\n");
      try
      {
        for (int i = 0; i < n; i++)
        {
          String s = "";
          for (int j = 0; j < n; j++)
          {
            if (j % m == 0) s += " |";
            for (int k = 0; k < n; k++)
            {
              double[] sol = x.Index(new int[] { i, j, k }).Level();
              if (sol[0] > 0.5)
              {
                s += " " + (k + 1);
                break;
              }
            }
          }
          Console.WriteLine(s + " |");
          if ((i + 1) % m == 0)
            Console.WriteLine("\n");
        }
      }
      catch (SolutionError se) { }
    }
  }
}

total_variation.cs

Listing 16.41 total_variation.cs Click here to download.
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File:      total_variation.cs
//
// Purpose:   Demonstrates how to solve a total
//            variation problem using the Fusion API.using System;
using System;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class total_variation
  {
    public static Model total_var(int n, int m) {
      Model M = new Model("TV");

      Variable u = M.Variable("u", new int[] {n + 1, m + 1}, Domain.InRange(0.0, 1.0));
      Variable t = M.Variable("t", new int[] {n, m}, Domain.Unbounded());

      // In this example we define sigma and the input image f as parameters
      // to demonstrate how to solve the same model with many data variants.
      // Of course they could simply be passed as ordinary arrays if that is not needed.
      Parameter sigma = M.Parameter("sigma");
      Parameter f = M.Parameter("f", n, m);

      Variable ucore = u.Slice(new int[] {0, 0}, new int[] {n, m});

      Expression deltax = Expr.Sub( u.Slice( new int[] {1, 0}, new int[] {n + 1, m} ), ucore );
      Expression deltay = Expr.Sub( u.Slice( new int[] {0, 1}, new int[] {n, m + 1} ), ucore );

      M.Constraint( Expr.Stack(2, t, deltax, deltay), Domain.InQCone().Axis(2) );

      M.Constraint( Expr.Vstack(sigma, Expr.Flatten( Expr.Sub(f,  ucore) ) ),
                    Domain.InQCone() );

      M.Objective( ObjectiveSense.Minimize, Expr.Sum(t) );

      return M;
    }

    public static void Main(string[] args)
    {
      Random randGen = new Random(0);
      
      int n = 100;
      int m = 200;
      double[] sigmas = { 0.0004, 0.0005, 0.0006 };

      // Create a parametrized model with given shape
      Model M = total_var(n, m);
      Parameter sigma = M.GetParameter("sigma");
      Parameter f     = M.GetParameter("f");
      Variable  ucore = M.GetVariable("u").Slice(new int[] {0,0}, new int[] {n,m});

      // Example: Linear signal with Gaussian noise    
      double[,] signal = new double[n,m];
      double[,] noise  = new double[n,m];
      double[,] fVal   = new double[n,m];
      double[,] sol    = new double[n,m];

      for(int i=0; i<n; i++) for(int j=0; j<m; j++) {
        signal[i,j] = 1.0*(i+j)/(n+m);
        noise[i,j] = Math.Sqrt(-2.0 * Math.Log(randGen.NextDouble())) * 
                     Math.Sin(2.0 * Math.PI * randGen.NextDouble()) *
                     0.08;
        fVal[i,j] = Math.Max( Math.Min(1.0, signal[i,j] + noise[i,j]), .0 );
      }
     
      // Set value for f
      f.SetValue(fVal);

      for(int iter=0; iter<3; iter++) {
        // Set new value for sigma and solve
        sigma.SetValue(sigmas[iter]*n*m);

        M.Solve();

        // Retrieve solution from ucore
        double[] ucoreLev = ucore.Level();
        for(int i=0; i<n; i++) for(int j=0; j<m; j++)
          sol[i,j] = ucoreLev[i*n+m];

        // Use the solution
        // ...

        Console.WriteLine("rel_sigma = " + sigmas[iter] + " total_var = " + M.PrimalObjValue());
      }

      M.Dispose();
    }
  }
}

tsp.cs

Listing 16.42 tsp.cs Click here to download.
//
//  Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
//  File:      tsp.cs
//
//  Purpose: Demonstrates a simple technique to the TSP
//           usign the Fusion API.
//
using System;
using System.Collections.Generic;
using System.Collections;
using mosek.fusion;

namespace mosek.fusion.example
{
  public class tsp {
    public static void Main(string [] args) {
      int []A_i = new int[] {0, 1, 2, 3, 1, 0, 2, 0};
      int []A_j = new int[] {1, 2, 3, 0, 0, 2, 1, 3};
      double []C_v =  new double[] {1.0, 1.0, 1.0, 1.0, 0.1, 0.1, 0.1, 0.1};

      int n = 4;
      Matrix adj = Matrix.Sparse(n, n, A_i, A_j, 1.0);
      Matrix costs = Matrix.Sparse(n, n, A_i, A_j, C_v);

      tsp_fusion(n, adj, costs , true, true);
      tsp_fusion(n, adj, costs , true, false);
    }

    public static void tsp_fusion(int n, Matrix A, Matrix C,
                                  bool remove_loops,
                                  bool remove_2_hop_loops) {
      Model M = new Model();

      Variable x = M.Variable(new int[] {n, n}, Domain.Binary());

      M.Constraint(Expr.Sum(x, 0), Domain.EqualsTo(1.0));
      M.Constraint(Expr.Sum(x, 1), Domain.EqualsTo(1.0));
      M.Constraint(x, Domain.LessThan( A ));

      M.Objective(ObjectiveSense.Minimize, Expr.Dot(C, x));

      if (remove_2_hop_loops)
        M.Constraint(Expr.Add(x, x.Transpose()), Domain.LessThan(1.0));

      if (remove_loops)
        M.Constraint(x.Diag(), Domain.EqualsTo(0.0));

      for (int it = 0; true; it++)
      {
        M.Solve();
        try
        {
          Console.WriteLine("\nit #{0} - solution cost: {1}", it, M.PrimalObjValue());
        } catch (SolutionError e) {}

        ArrayList cycles = new ArrayList();

        for (int i = 0; i < n; i++) {
          for (int j = 0; j < n; j++) {
            try {
              if ( x.Level()[i * n + j] <= 0.5) continue;
            } catch (SolutionError e) {}

            bool found = false;
            foreach (ArrayList c in cycles) {
              foreach (int[] a in c)
                if ( i == a[0] || j == a[0] ||
                     i == a[1] || j == a[1]  )
                {
                  c.Add( new int[] {i, j} );
                  found = true;
                  break;
                }
              if (found == true) break;
            }

            if (found != true) {
              ArrayList ll = new ArrayList();
              ll.Add( new int[] {i, j} );
              cycles.Add(ll);
            }
          }
        }

        Console.WriteLine("\ncycles:");
        foreach (ArrayList c in cycles) {
          foreach (int[] a in c)
            Console.Write("[{0},{1}] - ", a[0], a[1]);
          Console.Write("\n");
        }

        if (cycles.Count == 1) break;

        foreach (ArrayList c in cycles) {
          int [,] cc = new int[c.Count, 2];
          int i = 0;
          foreach (int[] a in c) {
            cc[i, 0] = a[0];
            cc[i, 1] = a[1];
            i++;
          }
          M.Constraint(Expr.Sum(x.Pick(cc)), Domain.LessThan( 1.0 * c.Count - 1 ));
        }
      }

      try {
        Console.WriteLine("\nsolution:");
        for (int i = 0; i < n; i++) {
          for (int j = 0; j < n; j++)
            Console.Write(" {0} ", ((int) x.Level()[i * n + j]));
          Console.Write("\n");
        }
      } catch (SolutionError e) {}
    }
  }
}