16 List of examples¶
List of examples shipped in the distribution of Fusion API for Java:
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.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: TrafficNetworkModel.java
//
// Purpose: Demonstrates a traffic 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.
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class TrafficNetworkModel {
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) throws SolutionError {
try(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;
int[] 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;
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.dot(e, 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();
return x.level();
}
}
public static void main(String[] args)
throws SolutionError {
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;
double[]flow = solve(n, source_idx, sink_idx,
arc_i, arc_j,
arc_sens,
arc_cap,
arc_base,
T);
System.out.println("Optimal flow:");
for (int i = 0; i < arc_i.length; ++i)
System.out.println("\tflow node" + arc_i[i] + "->node" + arc_j[i] + " = " + flow[arc_i[i] * n + arc_j[i]]);
}
}
alan.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: alan.java
//
// 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.
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class alan {
/////////////////////////////////////////////////////////////////////
// Problem data.
// Security names
private static final String[]
securities = { "hardware", "software", "show-biz", "t-bills" };
// Two sample mean returns on securities
private static final double[]
mean1 = { 9.0, 7.0, 11.0, 5.0 };
private static final double[]
mean2 = { 8.0, 9.0, 12.0, 7.0 };
// Target mean return
private static final double
target = 10.0;
private static final int numsec = securities.length;
// Factor of covariance matrix.
private static final double[][] U_data = {
{ 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 final Matrix U = Matrix.dense(U_data);
public static void main(String[] args)
throws SolutionError {
// Create a parametrized model
Model M = new Model("alan");
try {
Variable x = M.variable("x", numsec, Domain.greaterThan(0.0));
Variable t = M.variable("t", 1, Domain.greaterThan(0.0));
M.objective("minvar", ObjectiveSense.Minimize, t);
// 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("t > ||Ux||^2", Expr.vstack(0.5, t, Expr.mul(U, x)), Domain.inRotatedQCone());
// Solve with two values of mean
solve(M, mean, x, mean1);
solve(M, mean, x, mean2);
} finally {
M.dispose();
}
}
// Solve the model given some expected return vector
public static void solve(Model M, Parameter mean, Variable x, double[] meanVal)
throws SolutionError {
mean.setValue(meanVal);
System.out.print("Solve with mean = ");
for (int i = 0; i < numsec; ++i)
System.out.printf(" %e", meanVal[i]);
System.out.print("\n");
M.solve();
double[] solx = x.level();
System.out.print("Solution = ");
for (int i = 0; i < numsec; ++i)
System.out.printf(" %e", solx[i]);
System.out.print("\n");
}
}
baker.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: baker.java
//
// 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.
package com.mosek.fusion.examples;
import mosek.fusion.*;
import java.io.PrintWriter;
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)
throws SolutionError {
Matrix recipe = Matrix.dense(recipe_data);
Model M = new Model("Recipe");
try {
// "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(new PrintWriter(System.out));
// We solve and fetch the solution:
M.solve();
double[] res = production.level();
System.out.println("Solution:");
for (int i = 0; i < res.length; ++i) {
System.out.println(" Number of " + productnames[i] + " : " + res[i]);
}
System.out.println(" Revenue : $" +
(res[0] * revenue[0] + res[1] * revenue[1]));
} finally {
M.dispose();
}
}
}
breaksolver.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: breaksolver.java
//
// Purpose: Show how to break a long-running task.
//
// Requires a parameter defining a timeout in seconds.
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
import java.util.Random;
public class breaksolver {
private static void nshuffle(Random R, int[] a, int n) {
for (int i = 0; i < n; ++i) {
int idx = (int)(R.nextDouble() * (a.length - i)) + i;
int tmp = a[i]; a[i] = a[idx]; a[idx] = tmp;
}
}
public static void main(String[] args)
throws java.lang.InterruptedException, SolutionError {
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);
System.out.println("Build problem...");
try (Model M = new Model("SolveBinary")) {
M.setLogHandler(new java.io.PrintWriter(System.out));
Variable x = M.variable("x", n, Domain.binary());
M.objective(ObjectiveSense.Minimize, Expr.sum(x));
int[] idxs = new int[n];
for (int i = 0; i < n; ++i)
idxs[i] = i;
for (int i = 0; i < m; ++i) {
nshuffle(R, idxs, p);
M.constraint(Expr.sum(x.pick(java.util.Arrays.copyOf(idxs, p))),
Domain.equalsTo(p / 2));
}
System.out.println("Start thread...");
Thread T = new Thread() { public void run() { M.solve(); } };
T.start();
long T0 = System.currentTimeMillis();
while (true) {
if (System.currentTimeMillis() - T0 > timeout * 1000) {
System.out.println("Solver terminated due to timeout!");
M.breakSolver();
T.join();
break;
}
if (! T.isAlive()) {
System.out.println("Solver terminated before anything happened!");
T.join();
break;
}
}
} finally {
System.out.println("End.");
}
}
}
callback.java
/*
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
File : callback.java
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.
*/
package com.mosek.fusion.examples;
import mosek.*;
import mosek.fusion.*;
import java.util.Formatter;
import java.util.*;
import java.lang.*;
public class callback {
private static DataCallback makeUserCallback(final double maxtime,
final Model mod) {
return new DataCallback() {
public int callback(callbackcode caller,
double[] douinf,
int[] intinf,
long[] lintinf)
{
double opttime = 0.0;
int itrn;
double pobj, dobj, stime;
Formatter f = new Formatter(System.out);
switch (caller) {
case begin_intpnt:
f.format("Starting interior-point optimizer\n");
break;
case intpnt:
itrn = intinf[iinfitem.intpnt_iter.value ];
pobj = douinf[dinfitem.intpnt_primal_obj.value];
dobj = douinf[dinfitem.intpnt_dual_obj.value ];
stime = douinf[dinfitem.intpnt_time.value ];
opttime = douinf[dinfitem.optimizer_time.value ];
f.format("Iterations: %-3d\n", itrn);
f.format(" Elapsed time: %6.2f(%.2f)\n", opttime, stime);
f.format(" Primal obj.: %-18.6e Dual obj.: %-18.6e\n", pobj, dobj);
break;
case end_intpnt:
f.format("Interior-point optimizer finished.\n");
break;
case begin_primal_simplex:
f.format("Primal simplex optimizer started.\n");
break;
case update_primal_simplex:
itrn = intinf[iinfitem.sim_primal_iter.value ];
pobj = douinf[dinfitem.sim_obj.value ];
stime = douinf[dinfitem.sim_time.value ];
opttime = douinf[dinfitem.optimizer_time.value ];
f.format("Iterations: %-3d\n", itrn);
f.format(" Elapsed time: %6.2f(%.2f)\n", opttime, stime);
f.format(" Obj.: %-18.6e\n", pobj );
break;
case end_primal_simplex:
f.format("Primal simplex optimizer finished.\n");
break;
case begin_dual_simplex:
f.format("Dual simplex optimizer started.\n");
break;
case update_dual_simplex:
itrn = intinf[iinfitem.sim_dual_iter.value ];
pobj = douinf[dinfitem.sim_obj.value ];
stime = douinf[dinfitem.sim_time.value ];
opttime = douinf[dinfitem.optimizer_time.value ];
f.format("Iterations: %-3d\n", itrn);
f.format(" Elapsed time: %6.2f(%.2f)\n", opttime, stime);
f.format(" Obj.: %-18.6e\n", pobj);
break;
case end_dual_simplex:
f.format("Dual simplex optimizer finished.\n");
break;
case begin_bi:
f.format("Basis identification started.\n");
break;
case end_bi:
f.format("Basis identification finished.\n");
break;
default:
}
System.out.flush();
if (opttime >= maxtime) {
f.format("MOSEK is spending too much time. Terminate it.\n");
System.out.flush();
return 1;
}
return 0;
}
};
}
private static Progress makeProgressCallback() {
return new Progress() {
public int progress(callbackcode caller)
{
// Process the caller code here
return 0;
}
};
}
public static void main(String[] args) {
String slvr = "intpnt";
if (args.length < 1) {
System.out.println("Usage: ( psim | dsim | intpnt ) ");
}
if (args.length >= 1) slvr = args[0];
// Create a random 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();
try(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");
double maxtime = 0.07;
M.setDataCallbackHandler( makeUserCallback(maxtime, M) );
M.solve();
} catch (mosek.Exception e) {
System.out.println ("An error/warning was encountered");
System.out.println (e.toString());
throw e;
}
}
}
ceo1.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: ceo1.java
//
// 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)
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class ceo1 {
public static void main(String[] args)
throws SolutionError {
Model M = new Model("ceo1");
try {
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 conic exponential 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();
System.out.printf("x1,x2,x3 = %e, %e, %e\n", solx[0], solx[1], solx[2]);
// Get conic solution of expc
double[] expclvl = expc.level();
double[] expcsn = expc.dual();
System.out.printf("expc levels = %e", expclvl[0]);
for (int i = 1; i < expclvl.length; ++i)
System.out.printf(", %e", expclvl[i]);
System.out.print("\n");
System.out.printf("expc dual conic var levels = %e", expcsn[0]);
for (int i = 1; i < expcsn.length; ++i)
System.out.printf(", %e", expcsn[i]);
System.out.print("\n");
} finally {
M.dispose();
}
}
}
cqo1.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: cqo1.java
//
// 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 }
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class cqo1 {
public static void main(String[] args)
throws SolutionError {
Model M = new Model("cqo1");
try {
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, Domain.inQCone());
Constraint qc2 = M.constraint("qc2", z2, 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.writeTask("cqo1.task");
M.solve();
// Get the linear solution values
double[] solx = x.level();
double[] soly = y.level();
System.out.printf("x1,x2,x3 = %e, %e, %e\n", solx[0], solx[1], solx[2]);
System.out.printf("y1,y2,y3 = %e, %e, %e\n", soly[0], soly[1], soly[2]);
// Get conic solution of qc1
double[] qc1lvl = qc1.level();
double[] qc1sn = qc1.dual();
System.out.printf("qc1 levels = %e", qc1lvl[0]);
for (int i = 1; i < qc1lvl.length; ++i)
System.out.printf(", %e", qc1lvl[i]);
System.out.print("\n");
System.out.printf("qc1 dual conic var levels = %e", qc1sn[0]);
for (int i = 1; i < qc1sn.length; ++i)
System.out.printf(", %e", qc1sn[i]);
System.out.print("\n");
} finally {
M.dispose();
}
}
}
diet.java
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: diet.java
//
// Purpose: Solving Stigler's Nutrition model (DIET,SEQ=7)
//
// Source: GAMS Model library,
// Dantzig, G B, Chapter 27.1. In Linear Programming and Extensions.
// Princeton University Press, Princeton, New Jersey, 1963.
//
// Given a set of nutrients, foods with associated nutrient values, allowance of
// nutrients per day, the model will find the cheapest combination of foods
// which will provide the necessary nutrients.
//
// Arguments for construction:
// name - Model name.
// foods - List of M names of the foods.
// nutrients - List of N names of the nutrients.
// daily_allowance - List of N floats denoting the daily allowance of each
// nutrient.
// nutritive_value - Two-dimensional MxN array of floats where each row
// denotes the nutrient values for a single food per $ spent.
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class diet extends Model {
public class DataException extends RuntimeException {
public DataException(String msg) { super(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) {
super(name);
boolean 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));
dailyNutrients =
constraint("Nutrient Balance",
Expr.mul(nutrientValue, dailyPurchase),
Domain.greaterThan(dailyAllowance));
objective(ObjectiveSense.Minimize, Expr.sum(dailyPurchase));
finished = true;
} finally {
if (! finished)
dispose();
}
}
public double[] dailyPurchase()
throws SolutionError {
return dailyPurchase.level();
}
public double[] dailyNutrients()
throws SolutionError {
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, .6, 6, 55.8, .2, 0, 0, 0}, // margarine
{ 7.4, 448, 16.4, 19, 28.1, .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, .2, 0, .5, 5, 0}, // lard
{ 2.2, 333, .2, 139, 169.2, 6.4, 50.8, 316, 525}, // liver
{ 4.4, 249, .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., 70., .8, 12., 5., 1.8, 2.7, 18., 75. };
diet M = new diet("Stinglers Diet Problem",
foods,
nutrients,
daily_allowance,
nutritive_value);
try {
M.solve();
try {
System.out.println("Solution:");
double[] x = M.dailyPurchase();
for (int i = 0; i < x.length; ++i)
System.out.println(foods[i] + " : $" + x[i]);
System.out.println("Nutrients:");
double[] y = M.dailyNutrients();
for (int i = 0; i < y.length; ++i)
System.out.println(nutrients[i] + " : " + y[i] + " " + nutrient_unit[i] +
" (" + daily_allowance[i] + ")");
} catch (SolutionError e) {
System.out.println("Solution error : " + e.getMessage());
}
} finally {
M.dispose();
}
}
}
djc1.java
////
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: djc1.java
//
// 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
////
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class djc1 {
public static void main(String[] args)
throws SolutionError {
try(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(new java.io.PrintWriter(System.out)); // Enable log output
// Solve the problem
M.solve();
if (M.getPrimalSolutionStatus() == SolutionStatus.Optimal) {
double[] sol = x.level();
System.out.printf("[x0,x1,x2,x3] = [%e, %e, %e, %e]\n", sol[0], sol[1], sol[2], sol[3]);
}
else {
System.out.printf("Anoter solution status");
}
}
}
}
duality.java
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: duality.java
//
// Purpose: Show how to read the dual value of a constraint.
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class duality {
public static void main(String[] args)
throws SolutionError {
double[][] A = { { -0.5, 1.0 } };
double[] b = { 1.0 };
double[] c = { 1.0, 1.0 };
Model M = new Model("duality");
try {
Variable x = M.variable("x", 2, Domain.greaterThan(0.0));
Constraint con = M.constraint(Expr.sub(Expr.mul(Matrix.dense(A), x), b), Domain.equalsTo(0.0));
M.objective("obj", ObjectiveSense.Minimize, Expr.dot(c, x));
M.solve();
double[] xsol = x.level();
double[] ssol = x.dual();
double[] ysol = con.dual();
System.out.printf("x1,x2,s1,s2,y = %e, %e, %e, %e, %e\n", xsol[0], xsol[1], ssol[0], ssol[1], ysol[0]);
} finally {
M.dispose();
}
}
}
elastic.java
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: elastic.java
//
// Purpose: Demonstrates model parametrization on the example of an elastic net linear regression:
//
// min_x |Ax-b|_2 + lambda1*|x|_1 + lambda2*|x|_2
package com.mosek.fusion.examples;
import mosek.fusion.*;
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() throws SolutionError {
//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");
try {
// 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();
System.out.printf("Objective %.5f, solution %.3f, %.3f\n", M.primalObjValue(), sol[0], sol[1]);
// Increase lambda1
lambda1.setValue(0.5);
M.solve();
sol = x.level();
System.out.printf("Objective %.5f, solution %.3f, %.3f\n", 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();
System.out.printf("Objective %.5f, solution %.3f, %.3f\n", M.primalObjValue(), sol[0], sol[1]);
// And increase lamda2
lambda2.setValue(1.4145);
M.solve();
sol = x.level();
System.out.printf("Objective %.5f, solution %.3f, %.3f\n", M.primalObjValue(), sol[0], sol[1]);
}
catch (SolutionError e) {
System.out.println("Solution not available");
throw e;
}
finally {
M.dispose();
}
}
public static void main(String[] argv) throws SolutionError {
smallExample();
}
}
facility_location.java
////
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: facility_location.java
//
// 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
////
package com.mosek.fusion.examples;
import mosek.fusion.*;
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. } } );
static int N = customerloc.numRows();
public static void main(String[] args)
throws SolutionError {
Model M = new Model("FacilityLocation");
// Variable holding the facility location
try {
Variable f = M.variable("facility", new int[]{1, 2}, Domain.unbounded());
// Variable defining the euclidian distances to each customer
Variable d = M.variable("dist", new int[]{N, 1}, Domain.greaterThan(0.0));
// Variable defining the x and y differences to each customer;
Variable t = M.variable("t", new int[]{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();
double[] res = f.level();
System.out.println("Facility location = " + res[0] + "," + res[1]);
} finally {
M.dispose();
}
}
}
gp1.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: gp1.java
//
// 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
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
import java.lang.Math;
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.length;
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)
throws SolutionError
{
Model M = new Model("max_vol_box");
try {
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[] {Math.log(2.0/Aw), Math.log(2.0/Aw)});
M.constraint(Expr.dot(new double[] {0,1,1}, xyz), Domain.lessThan(Math.log(Af)));
M.constraint(Expr.dot(new double[] {1,-1,0}, xyz), Domain.inRange(Math.log(alpha),Math.log(beta)));
M.constraint(Expr.dot(new double[] {0,-1,1}, xyz), Domain.inRange(Math.log(gamma),Math.log(delta)));
M.setLogHandler(new java.io.PrintWriter(System.out));
M.solve();
double[] xyzVal = xyz.level();
double[] hwdVal = new double[3];
for(int i=0; i<3; i++) hwdVal[i] = Math.exp(xyzVal[i]);
return hwdVal;
} finally {
M.dispose();
}
}
public static void main(String[] args)
throws SolutionError
{
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);
System.out.format("h=%.4f w=%.4f d=%.4f\n", hwd[0], hwd[1], hwd[2]);
}
}
lo1.java
////
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: lo1.java
//
// 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
////
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class lo1 {
public static void main(String[] args)
throws SolutionError {
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'
try(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();
System.out.printf("[x0,x1,x2,x3] = [%e, %e, %e, %e]\n", sol[0], sol[1], sol[2], sol[3]);
}
}
}
logistic.java
////
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: logistic.java
//
// Purpose: Implements logistic regression with regulatization.
//
// Demonstrates using the exponential cone and log-sum-exp in Fusion.
package com.mosek.fusion.examples;
import mosek.fusion.*;
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,
boolean[] y,
double lamb)
{
int n = X.length;
int d = X[0].length; // 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)
throws SolutionError {
// Test: detect and approximate a circle using degree 2 polynomials
int n = 30;
double[][] X = new double[n*n][6];
boolean[] Y = new boolean[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.setLogHandler(new java.io.PrintWriter(System.out));
M.solve();
for(int i=0;i<6;i++)
System.out.println(theta.level()[i]);
}
}
lownerjohn_ellipsoid.java
package com.mosek.fusion.examples;
/**
Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
File: lownerjohn_ellipsoid.java
Purpose:
Computes the Lowner-John inner and outer ellipsoidal
approximations of a polytope.
References:
[1] "Lectures on Modern Optimization", Ben-Tal and Nemirovski, 2000.
[2] "MOSEK modeling manual", 2013
*/
import java.util.Arrays;
import mosek.fusion.*;
public class lownerjohn_ellipsoid {
/**
Purpose: Models the hypograph of the n-th power of the
determinant of a positive definite matrix. See [1,2] for more details.
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}
*/
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));
// (Z11*Z22*...*Znn) >= t^n
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;
}
/**
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.
*/
public static Object[] lownerjohn_inner(double[][] A, double[] b)
throws SolutionError {
Model M = new Model("lownerjohn_inner");
try {
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());
// (b-Ad, AC) generate cones
M.constraint("qc", Expr.hstack(Expr.sub(b, Expr.mul(A, d)), Expr.mul(A, C)),
Domain.inQCone());
// Objective: Maximize t
M.objective(ObjectiveSense.Maximize, t);
M.solve();
return new Object[] {C.level(), d.level()};
} finally {
M.dispose();
}
}
/**
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
maximize t
subject to t <= det(P)^(1/n)
|| P*xi - c ||_2 <= 1, i=1,...,m
P is PSD.
*/
public static Object[] lownerjohn_outer(double[][] x)
throws SolutionError {
Model M = new Model("lownerjohn_outer");
try {
int m = x.length;
int n = x[0].length;
//Setup variables
Variable t = M.variable("t", 1, Domain.greaterThan(0.0));
Variable P = det_rootn(M, t, n);
Variable c = M.variable("c", n, Domain.unbounded());
//(1, Px-c) in cone
M.constraint("qc",
Expr.hstack(Expr.ones(m),
Expr.sub(Expr.mul(x, P),
Var.reshape(Var.repeat(c, m), new int[] {m, n})
)
),
Domain.inQCone());
//Objective: Maximize t
M.objective(ObjectiveSense.Maximize, t);
M.solve();
return new Object[] {P.level(), c.level()};
} finally {
M.dispose();
}
}
/**************************************************************************************/
public static void main(String[] argv)
throws SolutionError {
//The list of vertices of a polygon
final double[][] p = { {0., 0.}, {1., 3.}, {5.5, 4.5}, {7., 4.}, {7., 1.}, {3., -2.} };
double[][] A = new double[p.length][];
double[] b = new double[p.length];
int n = p.length;
//The h-representation of that polygon
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];
}
Object[] ResInn = lownerjohn_inner(A, b);
Object[] ResOut = lownerjohn_outer(p);
double[] Ci = (double[]) ResInn[0];
double[] di = (double[]) ResInn[1];
double[] Po = (double[]) ResOut[0];
double[] co = (double[]) ResOut[1];
System.out.println("Inner:");
System.out.println(" C = " + Arrays.toString(Ci));
System.out.println(" d = " + Arrays.toString(di));
System.out.println("Outer:");
System.out.println(" P = " + Arrays.toString(Po));
System.out.println(" c = " + Arrays.toString(co));
}
}
lpt.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: lpt.java
//
// Purpose: Demonstrates how to solve the multi-processor
// scheduling problem using the Fusion API.
//
package com.mosek.fusion.examples;
import java.util.Random;
import java.util.*;
import java.io.*;
import java.lang.*;
import java.lang.Math.*;
import java.lang.Number.*;
import mosek.*;
import mosek.fusion.*;
public class lpt {
public static void main(String [] args) {
int n = 30; //Number of tasks
int m = 6; //Number of processors
double lb = 1.; //The range of lengths of short tasks
double ub = 5.;
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);
Arrays.sort(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.) );
M.constraint( Expr.sub( Var.repeat(t, m), Expr.mul(x, T)), Domain.greaterThan(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(new PrintWriter(System.out));
M.setSolverParam("mioTolRelGap", .01);
M.solve();
try {
System.out.printf("initial solution: \n");
for (int i = 0; i < m; i++) {
System.out.printf("M %d [", i);
for (int y = 0; y < n; y++)
System.out.printf( "%d, ", (int) init[i * n + y]);
System.out.printf("]\n");
}
System.out.print("MOSEK solution:\n");
for (int i = 0; i < m; i++) {
System.out.printf("M %d [", i);
for (int y = 0; y < n; y++)
System.out.printf( "%d, ", (int)x.index(i, y).level()[0]);
System.out.printf("]\n");
}
} catch (SolutionError e) {}
}
}
mico1.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: mico1.java
//
// 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
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class mico1 {
public static void main(String[] args)
throws SolutionError {
Model M = new Model("mico1");
try {
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();
System.out.printf("x, y = %f, %f\n", x.level()[0], y.level()[0]);
} finally {
M.dispose();
}
}
}
milo1.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: milo1.java
//
// Purpose: Demonstrates how to solve a small mixed
// integer linear optimization problem.
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class milo1 {
public static void main(String[] args)
throws SolutionError {
double[][] A = {
{ 50.0, 31.0 },
{ 3.0, -2.0 }
};
double[] c = { 1.0, 0.64 };
Model M = new Model("milo1");
try {
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();
System.out.printf("x1,x2 = %e, %e\n", sol[0], sol[1]);
System.out.printf("MIP rel gap = %.2f (%f)\n",
M.getSolverDoubleInfo("mioObjRelGap"),
M.getSolverDoubleInfo("mioObjAbsGap"));
} finally {
M.dispose();
}
}
}
mioinitsol.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: mioinitsol.java
//
// Purpose: Demonstrates how to solve a small mixed
// integer linear optimization problem.
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class mioinitsol {
public static void main(String[] args)
throws SolutionError {
int n = 4;
double[] c = { 7.0, 10.0, 1.0, 5.0 };
Model M = new Model("mioinitsol");
try {
Variable x = M.variable("x", n, Domain.greaterThan(0.0));
x.slice(0,3).makeInteger();
M.constraint(Expr.sum(x), Domain.lessThan(2.5));
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");
M.solve();
// Get the solution values
double[] sol = x.level();
System.out.printf("x = [");
for (int i = 0; i < n; i++) {
System.out.printf("%e, ", sol[i]);
}
System.out.printf("]\n");
// Was the initial solution used?
int constr = M.getSolverIntInfo("mioConstructSolution");
double constrVal = M.getSolverDoubleInfo("mioConstructSolutionObj");
System.out.println("Construct solution utilization: " + constr);
System.out.println("Construct solution objective: " + constrVal);
} finally {
M.dispose();
}
}
}
nearestcorr.java
/*
Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
File: nearestcorr.java
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
*/
package com.mosek.fusion.examples;
import mosek.*;
import mosek.fusion.*;
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];
int[] 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);
}
Matrix S = Matrix.sparse(N * (N + 1) / 2, N * N, msubi, msubj, mcof);
return Expr.mul(S, Expr.flatten(e));
}
/* Nearest correlation with nuclear norm penalty */
private static void nearestcorr_nn(Matrix A, double[] gammas, double[] res, int[] rank)
throws SolutionError {
int N = A.numRows();
Model M = new Model("NucNorm");
try {
// Setup variables
Variable t = M.variable("t", 1, Domain.unbounded());
Variable X = M.variable("X", Domain.inPSDCone(N));
Variable w = M.variable("w", N, Domain.greaterThan(0.0));
// (t, vec (X + diag(w) - A)) in Q
Expression 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() );
// Trace(X)
Expression TX = Expr.sum(X.diag());
for (int k = 0; k < gammas.length; ++k) {
// Objective: Minimize t + gamma*Tr(X)
M.objective(ObjectiveSense.Minimize, Expr.add(t, Expr.mul(gammas[k], TX)));
M.solve();
// Get the eigenvalues of X and approximate its rank
double[] d = new double[N];
LinAlg.syeig(mosek.uplo.lo, N, X.level(), d);
int rnk = 0; for (int i = 0; i < d.length; ++i) if (d[i] > 1e-6) ++rnk;
res[k] = t.level()[0];
rank[k] = rnk;
}
} finally {
M.dispose();
}
}
private static void nearestcorr(Matrix A)
throws SolutionError {
int N = A.numRows();
Model M = new Model("NearestCorrelation");
try {
// Setting up the variables
Variable X = M.variable("X", Domain.inPSDCone(N));
Variable 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
System.out.println("X = \n" + mattostr(X.level(), N));
System.out.println("t = \n" + mattostr(t.level(), N));
} finally {
M.dispose();
}
}
public static void main(String[] argv)
throws SolutionError {
int N = 5;
Matrix A = Matrix.dense(N, N,
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
});
double[] gammas = range(0.0, 2.0, 0.1);
double[] res = new double[gammas.length];
int[] rank = new int[gammas.length];
nearestcorr(A);
nearestcorr_nn(A, gammas, res, rank);
for (int i = 0; i < gammas.length; ++i)
System.out.printf("gamma = %.1f, res=%.3e, rank = %d\n", gammas[i], res[i], rank[i]);
}
/* Utility functions */
private static String mattostr(double[] a, int n) {
StringBuilder b = new StringBuilder();
java.util.Formatter f = new java.util.Formatter(b, java.util.Locale.US);
if ( a.length == 1) {
f.format("%f\n", a[0]);
return b.toString();
}
for (int i = 0; i < a.length; ++i) {
if ( i % n == 0 )
b.append("[");
f.format(" %.3e ", a[i]);
if ( ( i + 1 ) % n == 0 )
b.append("]\n");
}
return b.toString();
}
private 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.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: opt_server_sync.java
//
// Purpose : Demonstrates how to use MOSEK OptServer
// to solve optimization problem synchronously
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class opt_server_sync {
public static void main(String[] args) throws SolutionError {
if (args.length<1) {
System.out.println("Missing argument, syntax is:");
System.out.println(" java com.mosek.fusion.examples.opt_server_sync addr [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(new java.io.PrintWriter(System.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();
System.out.printf("x1,x2,x3 = %e, %e, %e\n", solx[0], solx[1], solx[2]);
}
}
parallel.java
/*
Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
File: parallel.java
Purpose: Demonstrates parallel optimization using solveBatch()
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;
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) throws SolutionError
{
// 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)
System.out.printf("Model %d: Status %s Solution Status %s objective %.3f time %.3f\n",
i,
status[i],
models[i].getPrimalSolutionStatus(),
models[i].primalObjValue(),
models[i].getSolverDoubleInfo("optimizerTime"));
else
System.out.printf("Model %d: not solved\n", i);
}
}
parameters.java
/*
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
File : parameters.java
Purpose : Demonstrates a very simple example about how to set
parameters and read information items
with MOSEK Fusion
*/
package com.mosek.fusion.examples;
import mosek.*;
import mosek.fusion.*;
public class parameters {
public static void main(String[] argv) {
// Create the Model
Model M = new Model();
System.out.println("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) {
System.out.println("Wrong parameter value");
}
// Define and solve an optimization problem here
// M.solve()
// After optimization:
System.out.println("Get MOSEK information items");
double tm = M.getSolverDoubleInfo("optimizerTime");
int it = M.getSolverIntInfo("intpntIter");
System.out.println("Time: " + tm);
System.out.println("Iterations: " + it);
}
}
pinfeas.java
////
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: pinfeas.java
//
// Purpose: Demonstrates how to fetch a primal infeasibility certificate
// for a linear problem
////
package com.mosek.fusion.examples;
import mosek.fusion.*;
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)
System.out.printf("%s[%d], dual = %e\n", name, i, duals[i]);
}
}
public static void main(String[] args) throws SolutionError {
// 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};
try(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(new java.io.PrintWriter(System.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
System.out.println("Variable bounds important for infeasibility: ");
analyzeCertificate("x", x.getSize(), x.dual(), eps);
// Go through constraint bounds
System.out.println("Constraint bounds important for infeasibility: ");
analyzeCertificate("s", s.getSize(), s.dual(), eps);
analyzeCertificate("d", d.getSize(), d.dual(), eps);
}
else {
System.out.println("The problem is not primal infeasible, no certificate to show");
}
}
}
}
portfolio_1_basic.java
/*
File : portfolio_1_basic.java
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
Purpose : Implements a basic portfolio optimization model.
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.ArrayList;
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)
throws mosek.fusion.SolutionError {
Model M = new Model("Basic Markowitz");
try {
// Redirect log output from the solver to stdout for debugging.
// if uncommented.
//M.setLogHandler(new java.io.PrintWriter(System.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/javafusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.format("Unexpected solution status: %s", solsta.toString()));
}
return dot(mu, x.level());
} finally {
M.dispose();
}
}
/*
The example. Reads in data and solves the portfolio models.
*/
public static void main(String[] argv)
throws java.io.IOException,
java.io.FileNotFoundException,
mosek.fusion.SolutionError {
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.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
{0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
{0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
{0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 },
{0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187},
{0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327},
{0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 }
};
System.out.println("\n-----------------------------------------------------------------------------------");
System.out.println("Basic Markowitz portfolio optimization");
System.out.println("-----------------------------------------------------------------------------------\n");
for ( int i = 0; i < gammas.length; ++i) {
double expret = BasicMarkowitz( n, mu, GT, x0, w, gammas[i]);
System.out.format("Expected return: %.4e Std. deviation: %.4e\n",
expret,
gammas[i]);
}
}
}
portfolio_2_frontier.java
/*
File : portfolio_2_frontier.java
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
Purpose : Implements a basic portfolio optimization model.
Computes points on the efficient frontier.
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.ArrayList;
import java.lang.Math;
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)
throws mosek.fusion.SolutionError {
Model M = new Model("Efficient frontier");
try {
//M.setLogHandler(new java.io.PrintWriter(System.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 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/javafusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.format("Unexpected solution status: %s", solsta.toString()));
}
frontier_mux[i] = dot(mu, x.level());
frontier_s[i] = s.level()[0];
}
} finally {
M.dispose();
}
}
/*
The example. Reads in data and solves the portfolio models.
*/
public static void main(String[] argv)
throws java.io.IOException,
java.io.FileNotFoundException,
mosek.fusion.SolutionError {
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.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
{0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
{0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
{0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 },
{0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187},
{0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327},
{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);
System.out.println("\n-----------------------------------------------------------------------------------");
System.out.println("Efficient frontier") ;
System.out.println("-----------------------------------------------------------------------------------\n");
System.out.format("%-12s %-12s %-12s\n", "alpha", "return", "std. dev.") ;
for (int i = 0; i < frontier_mux.length; ++i)
System.out.format("\t%-12.4f %-12.4e %-12.4e\n", alphas[i], frontier_mux[i], Math.sqrt(frontier_s[i]));
}
}
portfolio_3_impact.java
/*
File : portfolio_3_impact.java
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
Purpose : Implements a basic portfolio optimization model
with transaction costs of the form x^(3/2).
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.ArrayList;
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)
throws mosek.fusion.SolutionError {
Model M = new Model("Markowitz portfolio with market impact");
try {
//M.setLogHandler(new java.io.PrintWriter(System.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/javafusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.format("Unexpected solution status: %s", solsta.toString()));
}
if (xsol != null)
System.arraycopy(x.level(), 0, xsol, 0, n);
if (tsol != null)
System.arraycopy(t.level(), 0, tsol, 0, n);
} finally {
M.dispose();
}
}
/*
The example. Reads in data and solves the portfolio models.
*/
public static void main(String[] argv)
throws java.io.IOException,
java.io.FileNotFoundException,
mosek.fusion.SolutionError {
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);
System.out.println("\n-----------------------------------------------------------------------------------");
System.out.println("Markowitz portfolio optimization with market impact cost");
System.out.println("-----------------------------------------------------------------------------------\n");
System.out.format("Expected return: %.4e Std. deviation: %.4e Market impact cost: %.4e\n",
dot(mu, x),
gamma,
dot(m, t));
}
}
portfolio_4_transcost.java
/*
File : portfolio_4_transcost.java
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
Purpose : Implements a basic portfolio optimization model
with fixed setup costs and transaction costs
as a mixed-integer problem.
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.ArrayList;
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)
throws mosek.fusion.SolutionError {
// 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;
}
Model M = new Model("Markowitz portfolio with transaction costs");
try {
//M.setLogHandler(new java.io.PrintWriter(System.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/javafusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.format("Unexpected solution status: %s", solsta.toString()));
}
return x.level();
} finally {
M.dispose();
}
}
/*
The example. Reads in data and solves the portfolio models.
*/
public static void main(String[] argv)
throws java.io.IOException,
java.io.FileNotFoundException,
mosek.fusion.SolutionError {
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.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
{0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
{0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
{0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 },
{0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187},
{0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327},
{0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 }
};
double[] f = new double[n]; java.util.Arrays.fill(f, 0.01);
double[] g = new double[n]; java.util.Arrays.fill(g, 0.001);
double gamma = 0.36;
System.out.println("\n-----------------------------------------------------------------------------------");
System.out.println("Markowitz portfolio optimization with transaction cost");
System.out.println("-----------------------------------------------------------------------------------\n");
double[] x = new double[n];
x = MarkowitzWithTransactionsCost(n, mu, GT, x0, w, gamma, f, g);
System.out.println("Optimal portfolio: \n");
for ( int i = 0; i < x.length; ++i)
System.out.format("\tx[%-2d] %-12.4e\n", i, x[i]);
System.out.format("Expected return: %-12.4e", dot(mu, x));
}
}
portfolio_5_card.java
/*
File : portfolio_5_card.java
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
Description : Implements a basic portfolio optimization model
with cardinality constraints on number of assets traded.
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.ArrayList;
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: Maximum number of assets on 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)
throws mosek.fusion.SolutionError {
// 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;
}
Model M = new Model("Markowitz portfolio with cardinality bounds");
try {
//M.setLogHandler(new java.io.PrintWriter(System.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/javafusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.format("Unexpected solution status: %s", solsta.toString()));
}
double[] sol = x.level();
for(int j = 0; j < n; j++) results[i][j] = sol[j];
}
return results;
} finally {
M.dispose();
}
}
/*
The example. Reads in data and solves the portfolio models.
*/
public static void main(String[] argv)
throws java.io.IOException,
java.io.FileNotFoundException,
mosek.fusion.SolutionError {
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.25042, 0.09946, 0.09164, 0.06692, 0.08706, 0.09173, 0.08506},
{0. , 0. , 0.19914, 0.05867, 0.06453, 0.07367, 0.06468, 0.01914},
{0. , 0. , 0. , 0.20876, 0.04933, 0.03651, 0.09381, 0.07742},
{0. , 0. , 0. , 0. , 0.36096, 0.12574, 0.10157, 0.0571 },
{0. , 0. , 0. , 0. , 0. , 0.21552, 0.05663, 0.06187},
{0. , 0. , 0. , 0. , 0. , 0. , 0.22514, 0.03327},
{0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.2202 }
};
double gamma = 0.25;
int[] kValues = { 1, 2, 3, 4, 5, 6, 7, 8 };
double[][] results;
System.out.println("\n-----------------------------------------------------------------------------------");
System.out.println("Markowitz portfolio optimization with cardinality constraints");
System.out.println("-----------------------------------------------------------------------------------\n");
results = MarkowitzWithCardinality(n, mu, GT, x0, w, gamma, kValues);
for(int K = 1; K <= n; K++)
{
System.out.format("Bound %d Solution: ", K);
for (int i = 0; i < n; ++i)
System.out.format("%-12.4e ", results[K-1][i]);
System.out.println();
}
}
}
portfolio_6_factor.java
/*
File : portfolio_6_factor.java
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
Purpose : Implements a portfolio optimization model using factor model.
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;
import mosek.LinAlg;
import java.io.FileReader;
import java.io.BufferedReader;
import java.util.Arrays;
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.length;
int nj = m[0].length;
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.length;
int nj = m[0].length;
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.length;
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;
}
// Matrix multiplication
public static double[][] matrix_mul(double[][] a, double[][] b) {
int na = a.length;
int nb = b[0].length;
int k = b.length;
double[] vecm = new double[na * nb];
Arrays.fill(vecm, 0.0);
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 using factor model approximation of the covariance.
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)
throws mosek.fusion.SolutionError {
Model M = new Model("Factor model Markowitz");
try
{
// Redirect log output from the solver to stdout for debugging.
// if uncommented.
//M.setLogHandler(new java.io.PrintWriter(System.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/javafusion/accessing-solution.html about handling solution statuses.
throw new SolutionError(String.format("Unexpected solution status: %s", solsta.toString()));
}
return dot(mu, x.level());
}
finally
{
M.dispose();
}
}
/*
The example. Reads in data and solves the portfolio models.
*/
public static void main(String[] argv)
throws java.io.IOException,
java.io.FileNotFoundException,
mosek.fusion.SolutionError
{
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};
System.out.println("\n-----------------------------------------------------------------------------------");
System.out.println("Markowitz portfolio optimization based on a factor model.");
System.out.println("-----------------------------------------------------------------------------------\n");
for ( int i = 0; i < gammas.length; ++i)
{
double expret = FactorMarkowitz(n, mu, G_factor_T, theta, x0, w, gammas[i]);
System.out.format(
"Expected return: %.4e Std. deviation: %.4e\n",
expret,
gammas[i]
);
}
}
}
pow1.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: pow1.java
//
// 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
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class pow1 {
public static void main(String[] args)
throws SolutionError {
Model M = new Model("pow1");
try {
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 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();
System.out.printf("x, y, z = %e, %e, %e\n", solx[0], solx[1], solx[2]);
} finally {
M.dispose();
}
}
}
primal_svm.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: primal_svm.java
//
// Purpose: Implements a simple soft-margin SVM
// using the Fusion API.
package com.mosek.fusion.examples;
import java.util.Random;
import java.util.*;
import java.io.*;
import java.lang.*;
import java.lang.Math.*;
import java.lang.Number.*;
import mosek.fusion.*;
public class primal_svm {
private static final int nc = 10;
private static final int m = 50;
private static final int n = 3;
private static final int seed = 0;
private static final double var = 1.;
private static final double mean = 1.;
public static void main(String[] args)
throws SolutionError {
Random gen = new Random(seed);
int nump = gen.nextInt(m);
int numm = m - nump;
double [] y = new double[m];
Arrays.fill(y, 0, nump, 1.);
Arrays.fill(y, nump, m, -1.);
double [][] X = new double[m][n];
for (int i = 0; i < nump; i++)
for (int j = 0; j < n; j++)
X[i][j] = gen.nextGaussian() + 1.;
for (int i = nump; i < m; i++)
for (int j = 0; j < n; j++)
X[i][j] = gen.nextGaussian() - 1.;
Model M = new Model("Primal SVM");
try {
System.out.format("Number of data : %d\n", m);
System.out.format("Number of features: %d\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.));
M.constraint(
Expr.add( xi ,
Expr.mulElm( y,
Expr.sub( Expr.mul(X, w), Var.repeat(b, m) )
)
),
Domain.greaterThan( 1. )
);
M.constraint( Expr.vstack(1., t, w), Domain.inRotatedQCone());
System.out.println(" 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();
try {
System.out.format("%4f | %8f", c, b.level()[0] );
for (int j = 0; j < n; j++)
System.out.format(" | %8f", w.level()[j] );
System.out.print("\n");
} catch (FusionException e) {}
}
}
finally {
M.dispose();
}
}
}
qcqp_sdo_relaxation.java
/*
* Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
*
* File: qcqp_sdo_relaxation.java
*
* Purpose: Demonstrate how to use SDP to solve
* convex relaxation of a mixed-integer QCQP
**/
package com.mosek.fusion.examples;
import java.util.*;
import java.lang.*;
import java.lang.Math.*;
import java.io.*;
import mosek.*;
import mosek.fusion.*;
public class qcqp_sdo_relaxation {
public static void main(String [] args)
throws SolutionError {
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++)
A[i * n + j] = rnd.nextGaussian();
c[j] = rnd.nextDouble();
}
// P = A^T A
LinAlg.syrk(mosek.uplo.lo, mosek.transpose.yes,
n, m, 1.0, A, 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., q);
LinAlg.gemv(mosek.transpose.no, m, n, 1.0, A, c, 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] = java.lang.Math.round(xRound[i]);
double[] yRound = b.clone();
double[] xOpt = Mint.getVariable("x").level();
double[] yOpt = 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
System.out.println(M.getSolverDoubleInfo("optimizerTime") + " " + Mint.getSolverDoubleInfo("optimizerTime"));
System.out.println(java.lang.Math.sqrt(LinAlg.dot(m, yRound, yRound)) + " " +
java.lang.Math.sqrt(LinAlg.dot(m, yOpt, yOpt)));
}
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.) );
M.constraint( Z.index(n, n), Domain.equalsTo(1.) );
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|
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.java
/**
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
File : reoptimization.java
Purpose: Demonstrates how to solve a linear
optimization problem using the MOSEK API
and modify and re-optimize the problem.
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class reoptimization {
public static void printsol(double[] x) {
System.out.print("x = ");
for(double val : x)
System.out.printf("%.2f ", val);
System.out.println();
}
public static void main(String[] args)
throws SolutionError {
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
Model M = new Model();
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.java
/*
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
File : response.java
Purpose : This example demonstrates proper response handling
for problems solved with the interior-point optimizers.
*/
package com.mosek.fusion.examples;
import mosek.*;
import mosek.fusion.*;
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
/*
task.set_Stream(
streamtype.log,
new mosek.Stream()
{ public void stream(String msg) { System.out.print(msg); }});
*/
// (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();
System.out.print("Optimal value of x = ");
for(int i = 0; i < xsize; ++i)
System.out.print(xVal[i] + " ");
System.out.println("\nOptimal primal objective: " + M.primalObjValue());
// .. continue analyzing the solution
}
catch (OptimizeError e)
{
System.out.println("Optimization failed. Error: " + e.toString());
}
catch (SolutionError e)
{
// The solution with at least the expected status was not available.
// We try to diagnoze why.
System.out.println("Requested solution was not available.");
ProblemStatus prosta = M.getProblemStatus();
switch(prosta)
{
case DualInfeasible:
System.out.println("Dual infeasibility certificate found.");
break;
case PrimalInfeasible:
System.out.println("Primal infeasibility certificate found.");
break;
case Unknown:
// The solutions status is unknown. The termination code
// indicates why the optimizer terminated prematurely.
System.out.println("The solution status is unknown.");
StringBuffer symname = new StringBuffer();
StringBuffer desc = new StringBuffer();
Env.getcodedesc(rescode.fromValue(M.getSolverIntInfo("optimizeResponse")), symname, desc);
System.out.printf(" Termination code: %s %s\n", symname, desc);
break;
default:
System.out.println("Another unexpected problem status: " + prosta);
}
}
catch (java.lang.Exception e)
{
System.out.println("Unexpected error: " + e.toString());
}
M.dispose();
}
}
sdo1.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: sdo1.java
//
// 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)
//
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class sdo1 {
public static void main(String[] args) throws SolutionError {
Model M = new Model("sdo1");
try {
// 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();
System.out.println(java.util.Arrays.toString( X.level() ));
System.out.println(java.util.Arrays.toString( x.level() ));
} finally {
M.dispose();
}
}
}
sdo2.java
/*
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
File : sdo2.java
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.
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;
public class sdo2 {
public static void main(String[] args) throws SolutionError {
// 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);
Model M = new Model("sdo2");
try {
// 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(new java.io.PrintWriter(System.out));
M.solve();
// Print solution
System.out.println("Solution (vectorized):");
System.out.println(java.util.Arrays.toString( X1.level() ));
System.out.println(java.util.Arrays.toString( X2.level() ));
} finally {
M.dispose();
}
}
}
sdo3.java
/*
Copyright : Copyright (c) MOSEK ApS, Denmark. All rights reserved.
File : sdo3.java
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.
*/
package com.mosek.fusion.examples;
import mosek.fusion.*;
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) throws SolutionError {
// Sample input data
int n = 100;
int d = 4;
int k = 3;
double[] b = {9,10,11};
double[][][] A = new double[n*k][d][d];
for(int i=0; i<n*k; i++)
for(int s1=0; s1<d; s1++)
for(int s2=0; s2<=s1; s2++)
A[i][s1][s2] = A[i][s2][s1] = Math.random();
// Create a model with n semidefinite variables od dimension d x d
Model M = new Model("sdo3");
try {
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(new java.io.PrintWriter(System.out)); // Add logging
M.writeTask("sdo3.ptf"); // Save problem in readable format
M.solve();
// Get results. Each variable is a slice of X
System.out.println("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) {
System.out.println("X" + j + "=");
for(int s1=0; s1<d; s1++) {
for(int s2=0; s2<d; s2++)
System.out.print(Xj[s1*d+s1] + " ");
System.out.println();
}
}
}
}
finally {
M.dispose();
}
}
}
sospoly.java
////
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: sospoly.java
//
// 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.
package com.mosek.fusion.examples;
import mosek.fusion.*;
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), fill(i + 1, a));
else
return Matrix.sparse(n + 1, n + 1, range(n, i - n - 1, -1), range(i - n, n + 1), fill(2 * n + 1 - i, 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)
throws SolutionError {
Model M = new Model("smooth poly");
try {
int m = data.length;
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[data.length - 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))),
Domain.equalsTo(0.0));
nn_finite(M, lb, data[0][0], data[data.length - 1][0]);
M.objective(ObjectiveSense.Minimize, z);
M.solve();
return x.level();
} finally {
M.dispose();
}
}
public static void main(String[] argv)
throws SolutionError {
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);
System.out.println("fitpoly(data,2) -> " + arrtostr(x2));
System.out.println("fitpoly(data,4) -> " + arrtostr(x4));
System.out.println("fitpoly(data,8) -> " + 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) {
StringBuilder b = new StringBuilder("[");
java.util.Formatter f = new java.util.Formatter(b, java.util.Locale.US);
if (a.length > 0) f.format(", %.3e", a[0]);
for (int i = 1; i < a.length; ++i) f.format(", %.3e", a[i]);
b.append("]");
return b.toString();
}
}
sudoku.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: sudoku.java
//
// Purpose: A MILP-based SUDOKU solver
//
//
package com.mosek.fusion.examples;
import java.lang.*;
import java.util.*;
import java.io.*;
import mosek.fusion.*;
public class sudoku {
//fixed cells in human readable (i.e. 1-based) format
private static final 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}
};
public static void main(String[] args)
throws SolutionError {
int m = 3;
int n = m * m;
int nf = hr_fixed.length;
int [][] fixed = new int[nf][3];
for (int i = 0; i < nf; i++)
for (int d = 0; d < m; d++)
fixed[i][d] = hr_fixed[i][d] - 1;
try(Model M = new Model("SUDOKU")) {
M.setLogHandler(new PrintWriter(System.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.) );
//each number must appear 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++)
M.constraint( Expr.sum( x.slice( new int[] {i * m, j * m, k},
new int[] {(i + 1) * m, (j + 1) * m, k + 1}) ),
Domain.equalsTo(1.) );
M.constraint( x.pick(fixed) , Domain.equalsTo(1.0) ) ;
M.solve();
//print the solution, if any...
if ( M.getPrimalSolutionStatus() == SolutionStatus.Optimal )
print_solution(m, x);
else
System.out.println("No solution found!\n");
}
}
static void print_solution(int m, Variable x) {
int n = m * m;
System.out.println("\n");
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++)
try {
double sol[] = x.index( new int[] {i, j, k}).level();
if (sol[0] > 0.5) {
s += " " + (k + 1);
break;
}
} catch (SolutionError se) {}
}
System.out.println(s + " |");
if ((i + 1) % m == 0)
System.out.println("\n");
}
}
}
total_variation.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: total_variation.java
//
// Purpose: Demonstrates how to solve a total
// variation problem using the Fusion API.
//
package com.mosek.fusion.examples;
import java.util.Random;
import mosek.fusion.*;
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., 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)
throws SolutionError {
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] = randGen.nextGaussian() * 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
// ...
System.out.printf("rel_sigma = %f total_var = %.3f\n", sigmas[iter], M.primalObjValue());
}
M.dispose();
}
}
tsp.java
//
// Copyright: Copyright (c) MOSEK ApS, Denmark. All rights reserved.
//
// File: tsp.java
//
// Purpose: Demonstrates a simple technique to the TSP
// usign the Fusion API.
//
package com.mosek.fusion.examples;
import java.util.*;
import java.lang.*;
import java.io.*;
import java.lang.Math.*;
import mosek.*;
import mosek.fusion.*;
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., 1., 1., 1., 0.1, 0.1, 0.1, 0.1};
int n = 4;
Matrix costs = Matrix.sparse(n, n, A_i, A_j, C_v);
Matrix adj = Matrix.sparse(n, n, A_i, A_j, 1.);
tsp_fusion(n, adj, costs, true, false); //we do not remove loops of lenght 2
tsp_fusion(n, adj, costs, true, true); //we include all constraints
}
public static void tsp_fusion(int n,
Matrix A,
Matrix C,
boolean remove_loops,
boolean 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_loops)
M.constraint(x.diag(), Domain.equalsTo(0.));
if (remove_2_hop_loops)
M.constraint(Expr.add(x, x.transpose()), Domain.lessThan(1.0));
int it = 0;
while (true) {
it = it + 1;
M.solve();
try {
System.out.printf("\nit #%d - solution cost: %f\n", it, M.primalObjValue());
} catch (SolutionError e) {}
ArrayList<ArrayList< int[] > > cycles = new ArrayList<ArrayList< int[] > >();
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) {}
boolean found = false;
for (ArrayList<int[]> c : cycles) {
for (int[] a : 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< int[] > ll = new ArrayList<int[]>();
ll.add( new int[] {i, j} );
cycles.add(ll);
}
}
}
System.out.printf("\ncycles:\n");
for (ArrayList<int[]> c : cycles) {
for (int[] a : c)
System.out.printf("[%d,%d] - ", a[0], a[1]);
System.out.printf("\n");
}
if (cycles.size() == 1) break;
for (ArrayList<int []> c : cycles) {
int [][] cc = new int[c.size()][2];
c.toArray(cc);
M.constraint(Expr.sum(x.pick(cc)), Domain.lessThan( 1.0 * c.size() - 1 ));
}
}
try {
System.out.printf("\nsolution:\n");
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++)
System.out.printf(" %d ", ((int) x.level()[i * n + j]) );
System.out.printf("\n");
}
} catch (SolutionError e) {}
}
}
