16 List of examples¶
List of examples shipped in the distribution of Fusion API for .NET:
File |
Description |
|---|---|
Demonstrates a traffic network problem as a conic quadratic problem (CQO) |
|
A portfolio choice model |
|
A small bakery revenue maximization linear problem |
|
Shows how to break a long-running task |
|
An example of data/progress callback |
|
A simple conic exponential problem |
|
A simple conic quadratic problem |
|
Solving Stigler’s Nutrition model |
|
A simple problem with disjunctive constraints (DJC) |
|
Shows how to access the dual solution |
|
Linear regression with elastic net. Demonstrates model parametrization. |
|
Demonstrates a small one-facility location problem (CQO) |
|
A simple geometric program (GP) in conic form |
|
A simple linear problem |
|
Implements logistic regression and simple log-sum-exp (CEO) |
|
Computes the Lowner-John inner and outer ellipsoidal approximations of a polytope (SDO, Power Cone) |
|
Demonstrates how to solve the multi-processor scheduling problem and input an integer feasible point (MIP) |
|
A simple mixed-integer conic problem |
|
A simple mixed-integer linear problem |
|
A simple mixed-integer linear problem with an initial guess |
|
Solves the nearest correlation matrix problem (SDO, CQO) |
|
Uses MOSEK OptServer to solve an optimization problem synchronously |
|
Demonstrates parallel optimization using a batch method in MOSEK |
|
Shows how to set optimizer parameters and read information items |
|
Shows how to obtain and analyze a primal infeasibility certificate |
|
Portfolio optimization - basic Markowitz model |
|
Portfolio optimization - efficient frontier |
|
Portfolio optimization - market impact costs |
|
Portfolio optimization - transaction costs |
|
Portfolio optimization - cardinality constraints |
|
Portfolio optimization - factor model |
|
A simple power cone problem |
|
Implements a simple soft-margin Support Vector Machine (CQO) |
|
Demonstrate how to use SDP to solve convex relaxation of a mixed-integer QCQO problem |
|
Demonstrate how to modify and re-optimize a linear problem |
|
Demonstrates proper response handling |
|
A simple semidefinite problem with one matrix variable and a quadratic cone |
|
A simple semidefinite problem with two matrix variables |
|
A simple semidefinite problem with many matrix variables of the same dimension |
|
Models the cone of nonnegative polynomials and nonnegative trigonometric polynomials using Nesterov’s framework |
|
A SUDOKU solver (MIP) |
|
Demonstrates how to solve a total variation problem (CQO) |
|
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
//
// 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
//
// 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
// 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
//
// 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
//
// 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
//
// 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
//
// 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
/*
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
////
// 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
/*
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
// 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
////
// 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
//
// 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
////
// 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
////
// 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
/*
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
//
// 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 these lines to switch off feeding in the initial LPT solution
x.SetLevel(init);
M.SetSolverParam("mioConstructSol", "on");
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
//
// 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
//
// 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
//
// 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
/*
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, N, 1) );
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
//
// 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
/*
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
/*
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);
}
}
}
}
pinfeas.cs
////
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: pinfeas.cs
//
// Purpose: Demonstrates how to fetch a primal infeasibility certificate
// for a linear problem
////
using System;
using mosek.fusion;
namespace mosek.fusion.example
{
public class pinfeas {
//Analyzes and prints infeasibility certificate for a single object,
//which can be a variable or constraint
public static void analyzeCertificate(String name, // name of the analyzed object
long size, // size of the object
double[] duals, // actual dual values
double eps) { // tolerance determining when a dual value is considered important
for(int i = 0; i < size; i++) {
if (Math.Abs(duals[i]) > eps)
Console.WriteLine("{0}[{1}], dual = {2}", name, i, duals[i]);
}
}
public static void Main(string[] args) {
// Construct the sample model from the example in the manual
Matrix sMat = Matrix.Sparse(3, 7, new int[]{0,0,1,1,2,2,2},
new int[]{0,1,2,3,4,5,6},
new double[]{1,1,1,1,1,1,1});
double[] sBound = new double[]{200, 1000, 1000};
Matrix dMat = Matrix.Sparse(4, 7, new int[]{0,0,1,2,2,3,3},
new int[]{0,4,1,2,5,3,6},
new double[]{1,1,1,1,1,1,1});
double[] dBound = new double[]{1100, 200, 500, 500};
double[] c = new double[]{1, 2, 5, 2, 1, 2, 1};
using(Model M = new Model("pinfeas"))
{
Variable x = M.Variable("x", 7, Domain.GreaterThan(0));
Constraint s = M.Constraint("s", Expr.Mul(sMat, x), Domain.LessThan(sBound));
Constraint d = M.Constraint("d", Expr.Mul(dMat, x), Domain.EqualsTo(dBound));
M.Objective(ObjectiveSense.Minimize, Expr.Dot(c,x));
// Useful for debugging
M.WriteTask("pinfeas.ptf"); // Save to a readable file
M.SetLogHandler(Console.Out); // Enable log output
// Solve the problem
M.Solve();
// Check problem status
if (M.GetProblemStatus() == ProblemStatus.PrimalInfeasible) {
// Set the tolerance at which we consider a dual value as essential
double eps = 1e-7;
// We want to retrieve infeasibility certificates
M.AcceptedSolutionStatus(AccSolutionStatus.Certificate);
// Go through variable bounds
Console.WriteLine("Variable bounds important for infeasibility: ");
analyzeCertificate("x", x.GetSize(), x.Dual(), eps);
// Go through constraint bounds
Console.WriteLine("Constraint bounds important for infeasibility: ");
analyzeCertificate("s", s.GetSize(), s.Dual(), eps);
analyzeCertificate("d", d.GetSize(), d.Dual(), eps);
}
else {
Console.WriteLine("The problem is not primal infeasible, no certificate to show");
}
}
}
}
}
portfolio_1_basic.cs
/*
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();
// Check if the solution is an optimal point
SolutionStatus solsta = M.GetPrimalSolutionStatus();
if (solsta != SolutionStatus.Optimal)
{
// See https://docs.mosek.com/latest/dotnetfusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.Format("Unexpected solution status: {0}", solsta.ToString()));
}
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
/*
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();
// Check if the solution is an optimal point
SolutionStatus solsta = M.GetPrimalSolutionStatus();
if (solsta != SolutionStatus.Optimal)
{
// See https://docs.mosek.com/latest/dotnetfusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.Format("Unexpected solution status: {0}", solsta.ToString()));
}
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
/*
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();
// Check if the solution is an optimal point
SolutionStatus solsta = M.GetPrimalSolutionStatus();
if (solsta != SolutionStatus.Optimal)
{
// See https://docs.mosek.com/latest/dotnetfusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.Format("Unexpected solution status: {0}", solsta.ToString()));
}
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
/*
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();
// Check if the solution is an optimal point
SolutionStatus solsta = M.GetPrimalSolutionStatus();
if (solsta != SolutionStatus.Optimal)
{
// See https://docs.mosek.com/latest/dotnetfusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.Format("Unexpected solution status: {0}", solsta.ToString()));
}
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
/*
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();
// Check if the solution is an optimal point
SolutionStatus solsta = M.GetPrimalSolutionStatus();
if (solsta != SolutionStatus.Optimal)
{
// See https://docs.mosek.com/latest/dotnetfusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.Format("Unexpected solution status: {0}", solsta.ToString()));
}
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
/*
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;
}
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[] vector_sqrt(double[] m)
{
int ni = m.Length;
double[] sqrtm = new double[ni];
for (int i = 0; i < ni; ++i)
{
sqrtm[i] = Math.Sqrt(m[i]);
}
return sqrtm;
}
// 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
G_factor_T: The factor (dense) part of the factorized risk
theta: specific risk vector
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[,]G_factor_T,
double[] theta,
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(new Expression[] { Expr.ConstTerm(gamma),
Expr.Mul(G_factor_T, x),
Expr.MulElm(vector_sqrt(theta), x) }), Domain.InQCone());
// Solves the model.
M.Solve();
// Check if the solution is an optimal point
SolutionStatus solsta = M.GetPrimalSolutionStatus();
if (solsta != SolutionStatus.Optimal)
{
// See https://docs.mosek.com/latest/dotnetfusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.Format("Unexpected solution status: {0}", solsta.ToString()));
}
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[] theta = {0.0720, 0.0508, 0.0377, 0.0394, 0.0663, 0.0224, 0.0417, 0.0459};
double[,] P = cholesky(S_F);
double[,] G_factor = matrix_mul(B, P);
double[,] G_factor_T = transpose(G_factor);
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, G_factor_T, theta, x0, w, gamma);
Console.WriteLine("Expected return: {0,-12:f4} St deviation: {1,-12:f4} ", res, gamma);
}
}
}
}
pow1.cs
//
// 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
//
// 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
/*
* 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
/**
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
/*
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
//
// 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
/*
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
/*
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
////
// 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
//
// 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
//
// 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
//
// 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) {}
}
}
}
