diff --git a/Example2/n2-test.ipynb b/Example2/n2-test.ipynb index c4fb745..345b6cd 100644 --- a/Example2/n2-test.ipynb +++ b/Example2/n2-test.ipynb @@ -25,7 +25,7 @@ "np.set_printoptions(threshold=np.inf)\n", "\n", "# only show numpy output to five decimal places\n", - "np.set_printoptions(formatter={'float_kind':\"{:.5f}\".format})" + "np.set_printoptions(formatter={'float_kind':\"{:.8f}\".format})" ] }, { @@ -112,12 +112,12 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 1 rows, 2 columns and 2 nonzeros\n", - "Model fingerprint: 0x37f20dc3\n", + "Model fingerprint: 0xd1ce5e64\n", "Model has 2 quadratic objective terms\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", - " Objective range [5e-01, 1e+00]\n", - " QObjective range [1e+00, 1e+00]\n", + " Objective range [1e+00, 2e+00]\n", + " QObjective range [5e-01, 5e-01]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+00, 1e+00]\n", "Presolve time: 0.01s\n", @@ -133,19 +133,19 @@ "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 9.99499500e+05 -1.00100025e+06 2.00e+03 0.00e+00 1.00e+06 0s\n", - " 1 -4.99749639e-01 -9.99502094e+02 1.00e+00 3.41e-13 9.99e+02 0s\n", - " 2 -5.00116624e-01 -2.65997154e+02 1.00e-06 6.99e-15 1.33e+02 0s\n", - " 3 -5.00816469e-01 -1.01361177e+00 9.33e-10 0.00e+00 2.56e-01 0s\n", - " 4 -5.58837435e-01 -5.96696571e-01 3.43e-12 0.00e+00 1.89e-02 0s\n", - " 5 -5.62490657e-01 -5.64212774e-01 0.00e+00 0.00e+00 8.61e-04 0s\n", - " 6 -5.62500000e-01 -5.62501707e-01 0.00e+00 0.00e+00 8.53e-07 0s\n", - " 7 -5.62500000e-01 -5.62500002e-01 0.00e+00 0.00e+00 8.53e-10 0s\n", + " 0 -4.97498625e+05 5.00500125e+05 2.00e+03 0.00e+00 9.99e+05 0s\n", + " 1 2.66553322e+00 9.98756178e+02 1.11e+00 0.00e+00 1.05e+03 0s\n", + " 2 1.40062213e+00 2.93466464e+02 1.11e-06 0.00e+00 1.46e+02 0s\n", + " 3 1.40285882e+00 2.15287433e+00 1.35e-09 0.00e+00 3.75e-01 0s\n", + " 4 1.67958671e+00 1.77570144e+00 1.33e-15 0.00e+00 4.81e-02 0s\n", + " 5 1.74615363e+00 1.75032001e+00 0.00e+00 2.01e-16 2.08e-03 0s\n", + " 6 1.74999618e+00 1.75000030e+00 0.00e+00 0.00e+00 2.06e-06 0s\n", + " 7 1.75000000e+00 1.75000000e+00 0.00e+00 0.00e+00 2.06e-09 0s\n", "\n", "Barrier solved model in 7 iterations and 0.02 seconds (0.00 work units)\n", - "Optimal objective -5.62500000e-01\n", + "Optimal objective 1.75000000e+00\n", "\n", - "Standard: [0.25000 0.75000]\n" + "Standard: [0.00000001 0.99999999]\n" ] } ], @@ -164,8 +164,8 @@ "\n", "# define the objective functions for standard problem\n", "model_std.setObjective(\n", - " 0.5*w_std@(sigma@w_std) - lam*w_std.transpose()@mubar,\n", - "GRB.MINIMIZE)\n", + " w_std.transpose()@mubar - (lam/2)*w_std.transpose()@(sigma@w_std),\n", + "GRB.MAXIMIZE)\n", "\n", "# add sum-to-one constraint\n", "model_std.addConstr(np.ones([1,2]) @ w_std == 1, name=\"sum-to-one\")\n", @@ -196,44 +196,43 @@ "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", - "Optimize a model with 2 rows, 3 columns and 3 nonzeros\n", - "Model fingerprint: 0xb83c03fb\n", + "Optimize a model with 1 rows, 3 columns and 2 nonzeros\n", + "Model fingerprint: 0xea4b60b6\n", "Model has 2 quadratic objective terms\n", - "Model has 2 quadratic constraints\n", + "Model has 1 quadratic constraint\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " QMatrix range [1e-01, 4e+00]\n", - " Objective range [5e-01, 1e+00]\n", - " QObjective range [1e+00, 1e+00]\n", + " Objective range [1e+00, 2e+00]\n", + " QObjective range [5e-01, 5e-01]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+00, 1e+00]\n", - "Presolve removed 1 rows and 0 columns\n", "Presolve time: 0.01s\n", - "Presolved: 8 rows, 8 columns, 13 nonzeros\n", - "Presolved model has 3 second-order cone constraints\n", + "Presolved: 6 rows, 7 columns, 10 nonzeros\n", + "Presolved model has 2 second-order cone constraints\n", "Ordering time: 0.00s\n", "\n", "Barrier statistics:\n", - " AA' NZ : 1.900e+01\n", - " Factor NZ : 3.600e+01\n", - " Factor Ops : 2.040e+02 (less than 1 second per iteration)\n", + " AA' NZ : 1.500e+01\n", + " Factor NZ : 2.100e+01\n", + " Factor Ops : 9.100e+01 (less than 1 second per iteration)\n", " Threads : 1\n", "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 -6.30136986e-01 -6.30136986e-01 1.61e+00 8.66e-01 7.30e-01 0s\n", - " 1 -1.14616869e+00 -1.60124001e+00 4.19e-01 9.52e-07 1.64e-01 0s\n", - " 2 -6.51626451e-01 -8.34836457e-01 4.61e-07 1.05e-12 1.67e-02 0s\n", - " 3 -6.72758817e-01 -6.83613315e-01 2.69e-08 3.87e-14 9.87e-04 0s\n", - " 4 -6.73595613e-01 -6.73729333e-01 2.55e-14 2.00e-15 1.22e-05 0s\n", - " 5 -6.73606055e-01 -6.73616513e-01 1.05e-12 9.30e-13 9.51e-07 0s\n", - " 6 -6.73611064e-01 -6.73611270e-01 6.32e-11 1.09e-11 1.88e-08 0s\n", + " 0 -5.73852130e-02 1.22656250e+00 8.97e-01 1.41e-01 3.51e-01 0s\n", + " 1 5.91715790e-01 9.52821847e-01 2.13e-01 1.51e-07 6.38e-02 0s\n", + " 2 5.15372316e-01 5.96285470e-01 2.34e-07 6.30e-10 8.99e-03 0s\n", + " 3 5.49225634e-01 5.53202935e-01 4.81e-10 1.05e-15 4.42e-04 0s\n", + " 4 5.52127256e-01 5.52382277e-01 1.12e-13 5.55e-16 2.83e-05 0s\n", + " 5 5.52206847e-01 5.52232334e-01 3.52e-12 3.66e-15 2.83e-06 0s\n", + " 6 5.52216328e-01 5.52218236e-01 1.38e-11 2.96e-14 2.12e-07 0s\n", "\n", "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n", - "Optimal objective -6.73611064e-01\n", + "Optimal objective 5.52216328e-01\n", "\n", - "Standard: [0.25000 0.75000]\n", - "Robust: [0.41660 0.58340]\n" + "Standard: [0.00000001 0.99999999]\n", + "Robust: [0.83155134 0.16844866]\n" ] } ], @@ -245,18 +244,17 @@ "model_rbs = gp.Model(\"n02robust\")\n", "\n", "# initialise w for both models, z for robust model \n", - "w_rbs = model_rbs.addMVar(shape=n, vtype=GRB.CONTINUOUS, name=\"w\")\n", - "z_rbs = model_rbs.addVar(name=\"z\")\n", + "w_rbs = model_rbs.addMVar(shape=n, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "z_rbs = model_rbs.addVar(lb=0.0, name=\"z\")\n", "\n", "model_rbs.setObjective(\n", - " 0.5*w_rbs@(sigma@w_rbs) - lam*w_rbs.transpose()@mubar - kappa*z_rbs,\n", - "GRB.MINIMIZE)\n", + " w_rbs.transpose()@mubar - (lam/2)*w_rbs.transpose()@(sigma@w_rbs) - kappa*z_rbs,\n", + "GRB.MAXIMIZE)\n", "\n", "# add sum-to-half constraints to both models\n", "model_rbs.addConstr(np.ones([1,2]) @ w_std == 1, name=\"sum-to-one\")\n", "# add quadratic uncertainty constraint to the robust model\n", - "model_rbs.addConstr(z_rbs**2 <= np.inner(w_rbs, omega@w_rbs), name=\"uncertainty\")\n", - "model_rbs.addConstr(z_rbs >= 0, name=\"z positive\")\n", + "model_rbs.addConstr(z_rbs**2 >= w_rbs.transpose()@omega@w_rbs, name=\"uncertainty\")\n", "\n", "# solve problem with Gurobi, print alongside standard solution for comparison\n", "model_rbs.optimize()\n", @@ -268,7 +266,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that if $\\kappa = 0$, the standard solution not _exactly_ recovered, only accurate to four decimal places." + "Note that if $\\kappa = 0$, the standard solution not _exactly_ recovered, only accurate to seven decimal places." ] }, { @@ -285,44 +283,42 @@ "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", - "Optimize a model with 2 rows, 3 columns and 3 nonzeros\n", - "Model fingerprint: 0x48b37c4b\n", + "Optimize a model with 1 rows, 3 columns and 2 nonzeros\n", + "Model fingerprint: 0xcd1e4542\n", "Model has 2 quadratic objective terms\n", - "Model has 2 quadratic constraints\n", + "Model has 1 quadratic constraint\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " QMatrix range [1e-01, 4e+00]\n", - " Objective range [5e-01, 1e+00]\n", - " QObjective range [1e+00, 1e+00]\n", + " Objective range [1e+00, 2e+00]\n", + " QObjective range [5e-01, 5e-01]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+00, 1e+00]\n", - "Presolve removed 1 rows and 0 columns\n", "Presolve time: 0.01s\n", - "Presolved: 8 rows, 8 columns, 13 nonzeros\n", - "Presolved model has 3 second-order cone constraints\n", + "Presolved: 6 rows, 7 columns, 10 nonzeros\n", + "Presolved model has 2 second-order cone constraints\n", "Ordering time: 0.00s\n", "\n", "Barrier statistics:\n", - " AA' NZ : 1.900e+01\n", - " Factor NZ : 3.600e+01\n", - " Factor Ops : 2.040e+02 (less than 1 second per iteration)\n", + " AA' NZ : 1.500e+01\n", + " Factor NZ : 2.100e+01\n", + " Factor Ops : 9.100e+01 (less than 1 second per iteration)\n", " Threads : 1\n", "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 -6.30136986e-01 -6.30136986e-01 1.67e+00 4.15e-01 3.76e-01 0s\n", - " 1 -7.63463502e-01 -9.90451466e-01 3.34e-01 4.57e-07 6.75e-02 0s\n", - " 2 -5.29769678e-01 -6.77289143e-01 3.67e-07 5.02e-13 1.34e-02 0s\n", - " 3 -5.47218963e-01 -5.75552518e-01 4.04e-13 3.33e-16 2.58e-03 0s\n", - " 4 -5.61928380e-01 -5.63204951e-01 4.66e-15 1.11e-16 1.16e-04 0s\n", - " 5 -5.62460343e-01 -5.62501958e-01 9.44e-14 2.22e-16 3.78e-06 0s\n", - " 6 -5.62499955e-01 -5.62500028e-01 2.26e-11 6.10e-14 6.60e-09 0s\n", + " 0 1.22656250e+00 1.22656250e+00 1.08e+00 6.94e-01 4.95e-01 0s\n", + " 1 1.61470687e+00 1.61496484e+00 5.67e-02 1.55e-01 6.97e-02 0s\n", + " 2 1.67968985e+00 1.70168242e+00 6.24e-08 2.02e-02 1.10e-02 0s\n", + " 3 1.74763426e+00 1.75020504e+00 7.06e-14 2.52e-05 2.97e-04 0s\n", + " 4 1.74999548e+00 1.75001176e+00 1.87e-14 2.84e-08 1.82e-06 0s\n", + " 5 1.74999978e+00 1.75000067e+00 4.92e-13 9.33e-10 9.90e-08 0s\n", "\n", - "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n", - "Optimal objective -5.62499955e-01\n", + "Barrier solved model in 5 iterations and 0.03 seconds (0.00 work units)\n", + "Optimal objective 1.74999978e+00\n", "\n", - "Standard: [0.2500000031 0.7499999969]\n", - "Robust: [0.2499847190 0.7500152810]\n" + "Standard: [0.0000000076 0.9999999924]\n", + "Robust: [0.0000001531 0.9999998469]\n" ] } ], @@ -334,18 +330,17 @@ "model_rbs = gp.Model(\"n02robust_k0\")\n", "\n", "# initialise w for both models, z for robust model \n", - "w_rbs = model_rbs.addMVar(shape=n, vtype=GRB.CONTINUOUS, name=\"w\")\n", - "z_rbs = model_rbs.addVar(name=\"z\")\n", + "w_rbs = model_rbs.addMVar(shape=n, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "z_rbs = model_rbs.addVar(lb=0.0, name=\"z\")\n", "\n", "model_rbs.setObjective(\n", - " 0.5*w_rbs@(sigma@w_rbs) - lam*w_rbs.transpose()@mubar - kappa*z_rbs,\n", - "GRB.MINIMIZE)\n", + " w_rbs.transpose()@mubar - (lam/2)*w_rbs.transpose()@(sigma@w_rbs) - kappa*z_rbs,\n", + "GRB.MAXIMIZE)\n", "\n", "# add sum-to-half constraints to both models\n", "model_rbs.addConstr(np.ones([1,2]) @ w_std == 1, name=\"sum-to-one\")\n", "# add quadratic uncertainty constraint to the robust model\n", - "model_rbs.addConstr(z_rbs**2 <= np.inner(w_rbs, omega@w_rbs), name=\"uncertainty\")\n", - "model_rbs.addConstr(z_rbs >= 0, name=\"z positive\")\n", + "model_rbs.addConstr(z_rbs**2 >= w_rbs.transpose()@omega@w_rbs, name=\"uncertainty\")\n", "\n", "# solve problem with Gurobi, print alongside standard solution for comparison\n", "model_rbs.optimize()\n", diff --git a/Scratch/numerical-error.ipynb b/Scratch/numerical-error.ipynb index 63583a3..56a658e 100644 --- a/Scratch/numerical-error.ipynb +++ b/Scratch/numerical-error.ipynb @@ -6,7 +6,16 @@ "source": [ "# Numerical Error\n", "\n", - "The purpose of this notebook is to illustrate the numerical error which comes from setting $\\kappa = 0$." + "The purpose of this notebook is to illustrate the numerical error which comes from setting $\\kappa = 0$ and using \n", + "\n", + "```python3\n", + "model.addConstr(z**2 <= np.inner(w, omega@w), name=\"uncertainty\")\n", + "```\n", + "as the cone constraint rather than\n", + "```python3\n", + "model.addConstr(z**2 <= w.transpose()@omega@w, name=\"uncertainty\")\n", + "```\n", + "This has been fixed in the latest version of all other notebooks." ] }, { @@ -162,7 +171,7 @@ "model = gp.Model(\"standardGS\")\n", "\n", "# define variable of interest as a continuous \n", - "w = model.addMVar(shape=problem_size, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "w = model.addMVar(shape=problem_size, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\")\n", "\n", "# set the objective function\n", "model.setObjective(\n", @@ -172,8 +181,8 @@ "# add sub-to-half constraints\n", "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", "# add weight-bound constraints\n", - "model.addConstr(w >= lower_bound, name=\"lower bound\")\n", - "model.addConstr(w <= upper_bound, name=\"upper bound\")\n", + "model.addConstr(w >= lower_bound, name=\"lower-bound\")\n", + "model.addConstr(w <= upper_bound, name=\"upper-bound\")\n", "\n", "# solve the problem with Gurobi\n", "model.optimize()\n", @@ -231,8 +240,8 @@ "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", - "Optimize a model with 23 rows, 4 columns and 25 nonzeros\n", - "Model fingerprint: 0x97227d81\n", + "Optimize a model with 22 rows, 4 columns and 24 nonzeros\n", + "Model fingerprint: 0x52c6cd50\n", "Model has 3 quadratic objective terms\n", "Model has 3 quadratic constraints\n", "Coefficient statistics:\n", @@ -242,7 +251,7 @@ " QObjective range [1e+00, 5e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [5e-01, 1e+00]\n", - "Presolve removed 22 rows and 1 columns\n", + "Presolve removed 21 rows and 1 columns\n", "Presolve time: 0.01s\n", "Presolved: 9 rows, 8 columns, 14 nonzeros\n", "Presolved model has 3 second-order cone constraints\n", @@ -263,7 +272,7 @@ " 4 -8.12263522e-01 -8.12579837e-01 1.20e-14 2.22e-16 2.64e-05 0s\n", " 5 -8.12499544e-01 -8.12500834e-01 8.52e-13 3.00e-15 1.07e-07 0s\n", "\n", - "Barrier solved model in 5 iterations and 0.02 seconds (0.00 work units)\n", + "Barrier solved model in 5 iterations and 0.03 seconds (0.00 work units)\n", "Optimal objective -8.12499544e-01\n", "\n", "Standard: w = [0.50000 0.37500 0.12500],\n", @@ -277,8 +286,8 @@ "model = gp.Model(\"robustGS\")\n", "\n", "# define variables of interest as a continuous\n", - "w = model.addMVar(shape=problem_size, vtype=GRB.CONTINUOUS, name=\"w\")\n", - "z = model.addVar(name=\"z\")\n", + "w = model.addMVar(shape=problem_size, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "z = model.addVar(lb=0.0, name=\"z\")\n", "\n", "# setup the robust objective function\n", "model.setObjective(\n", @@ -287,12 +296,12 @@ "\n", "# add quadratic uncertainty constraint\n", "model.addConstr(z**2 <= np.inner(w, omega@w), name=\"uncertainty\")\n", - "model.addConstr(z >= 0, name=\"z positive\")\n", + "# model.addConstr(z**2 <= w.transpose()@omega@w, name=\"uncertainty\")\n", "# add sub-to-half constraints\n", "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", "# add weight-bound constraints~\n", - "model.addConstr(w >= lower_bound, name=\"lower bound\")\n", - "model.addConstr(w <= upper_bound, name=\"upper bound\")\n", + "model.addConstr(w >= lower_bound, name=\"lower-bound\")\n", + "model.addConstr(w <= upper_bound, name=\"upper-bound\")\n", "\n", "# solve the problem with Gurobi\n", "model.optimize()\n", @@ -303,7 +312,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that even though $\\kappa = 0$, we have $z\\neq0$. This suggests that a small amount of numerical error has been introduced as a result of asking Gurobi to optimize $z$ as well as $w$. " + "Note that even though $\\kappa = 0$, we have $z\\neq0$. This suggests that a small amount of numerical error has been introduced." ] } ], diff --git a/Scratch/risk-minimization.ipynb b/Scratch/risk-minimization.ipynb new file mode 100644 index 0000000..7a6d8c0 --- /dev/null +++ b/Scratch/risk-minimization.ipynb @@ -0,0 +1,824 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Alternative Framing\n", + "\n", + "In the [main notebook](../robust-genetics.ipynb) it's mentioned that the mathematics is more complicated when the robust optimization problem is framed as a risk minimization problem, rather than return maximization as done there. This notebook includes the calculations of the alternative framing to illustrate that point." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np # defines matrix structures\n", + "from qpsolvers import solve_qp # used for quadratic optimization\n", + "import gurobipy as gp # Gurobi optimization interface (1)\n", + "from gurobipy import GRB # Gurobi optimization interface (2)\n", + "\n", + "# want to round rather than truncate when printing\n", + "np.set_printoptions(threshold=np.inf)\n", + "\n", + "# only show numpy output to five decimal places\n", + "np.set_printoptions(formatter={'float_kind':\"{:.5f}\".format})" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Standard Problem\n", + "\n", + "In the context of genetic selection, we want to maximize selection of genetic merit (a measure of desirable traits) while minimizing risks due to inbreeding. This can be formed mathematically as a problem about minimizing risk for varying levels of return. We say\n", + "$$\n", + " \\min_w \\left(\\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu\\right)\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u,\n", + "$$\n", + "where $w$ is the vector of proportional contributions, $\\Sigma$ is a matrix encoding risk, $\\mu$ is a vector encoding returns, $l$ encodes lower bounds on contributions, $u$ encodes upper bounds on contributions, and $M$ is matrix which encodes whether each candidate is a sire or dam as described previously.\n", + "\n", + "In this representation of the problem, $\\lambda$ is a control variable which balances how we trade of between risk and return. Each value of $\\lambda$ will give a different solution on the critical frontier of the problem." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Toy Example ($n = 3$)\n", + "\n", + "Lets see how this works in an example. We will start by looking how this problem might be solving using Python's [qpsolvers](https://qpsolvers.github.io/qpsolvers/index.html) library. Consider the problem where\n", + "$$\n", + " \\mu = \\begin{bmatrix} 1 \\\\ 5 \\\\ 2 \\end{bmatrix},\\quad\n", + " \\Sigma = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 5 & 0 \\\\ 0 & 0 & 3 \\end{bmatrix}, \\quad\n", + " \\mathcal{S} = \\lbrace 1 \\rbrace, \\quad\n", + " \\mathcal{D} = \\lbrace 2, 3 \\rbrace, \\quad\n", + " l = {\\bf 0}, \\quad\n", + " u = {\\bf 1}.\n", + "$$\n", + "We define these variables in Python using the following code." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# KEY PROBLEM VARIABLES\n", + "problem_size = 3\n", + "expected_breeding_values = np.array([\n", + " 1.0,\n", + " 5.0,\n", + " 2.0\n", + "])\n", + "relationship_matrix = np.array([\n", + " [1, 0, 0],\n", + " [0, 5, 0],\n", + " [0, 0, 3]\n", + "])\n", + "sire_indices = [0]\n", + "dam_indices = [1,2]\n", + "lower_bound = np.full((problem_size, 1), 0.0)\n", + "upper_bound = np.full((problem_size, 1), 1.0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We have the additional variables which need setting up so that the problem works in `qpsolvers`. " + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "# OPTIMIZATION SETUP VARIABLES\n", + "lam = 0.5\n", + "# define the M so that column i is [1;0] if i is a sire and [0;1] otherwise \n", + "M = np.zeros((2, problem_size))\n", + "M[0, sire_indices] = 1\n", + "M[1, dam_indices] = 1\n", + "# define the right hand side of the constraint Mx = m\n", + "m = np.array([[0.5], [0.5]])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Finally, we solve the problem using the modules' `solve_qp` function. This utilises Gurobi via an API, a fact which will be important once we start to consider larger problem sizes. " + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Set parameter Username\n", + "Academic license - for non-commercial use only - expires 2025-02-26\n", + "QP solution: w = [0.50000 0.37500 0.12500]\n" + ] + } + ], + "source": [ + "# SOLVE THE PROBLEM\n", + "def solveGS(lam):\n", + " return solve_qp(\n", + " P = relationship_matrix,\n", + " q = -lam*expected_breeding_values,\n", + " G = None,\n", + " h = None,\n", + " A = M,\n", + " b = m,\n", + " lb = lower_bound,\n", + " ub = upper_bound,\n", + " solver = \"gurobi\"\n", + " )\n", + "\n", + "print(f\"QP solution: w = {solveGS(lam)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Excellent, we have code which is telling us that our optimal contributions for $\\lambda = 0.5$ are $w_1 = 0.5$, $w_2 = 0.375$, and $w_3 = 0.125$. We could vary our $\\lambda$ value too and find other points on the frontier." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "lambda = 0.0: w = [0.50000 0.18750 0.31250]\n", + "lambda = 0.2: w = [0.50000 0.26250 0.23750]\n", + "lambda = 0.4: w = [0.50000 0.33750 0.16250]\n", + "lambda = 0.6: w = [0.50000 0.41250 0.08750]\n", + "lambda = 0.8: w = [0.50000 0.48750 0.01250]\n", + "lambda = 1.0: w = [0.50000 0.50000 0.00000]\n" + ] + } + ], + "source": [ + "print(f\"lambda = 0.0: w = {solveGS(0.0)}\")\n", + "print(f\"lambda = 0.2: w = {solveGS(0.2)}\")\n", + "print(f\"lambda = 0.4: w = {solveGS(0.4)}\")\n", + "print(f\"lambda = 0.6: w = {solveGS(0.6)}\")\n", + "print(f\"lambda = 0.8: w = {solveGS(0.8)}\")\n", + "print(f\"lambda = 1.0: w = {solveGS(1.0)}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Robust Optimization\n", + "\n", + "As discussed, rather than treating $\\mu$ as known we must model it using a univariate normal distribution, $\\mu\\sim N(\\bar{\\mu}, \\Omega)$. We adjust the objective function to model the inherent uncertainty in the problem, and in the risk-minimization framing the problem the bilevel optimization problem is\n", + "$$\n", + " \\min_w \\left(\\frac{1}{2}w^{T}\\Sigma w - \\lambda\\max_{\\mu\\in U_{\\mu}} \\left(w^{T}\\mu \\right)\\right)\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + "$$\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Quadratic Uncertainty Sets\n", + "\n", + "We define $U_\\mu$ as a quadratic uncertainty set,\n", + "$$\n", + " U_{\\mu} := \\left\\lbrace \\mu :\\ {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\leq \\kappa^2 \\right\\rbrace.\n", + "$$\n", + "This means that our lower level problem $\\max_{\\mu\\in U_{\\mu}} w^{T}\\mu$ becomes\n", + "$$\n", + " \\max_{\\mu} w^{T}\\mu\\quad \\text{subject to}\\quad {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\leq \\kappa^2,\n", + "$$\n", + "or, in standard form,\n", + "$$\n", + " \\min_{\\mu} \\left(-w^{T}\\mu\\right) \\quad\\text{subject to}\\quad \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\geq 0.\n", + "$$\n", + "This is convex, since $\\Omega$ is a positive definite matrix and $-w^{T}\\mu$ is an affine function. If we define a Lagrangian multiplier $\\rho\\in\\mathbb{R}$, the KKT conditions of this problem are:\n", + "\\begin{align}\n", + " \\nabla_{\\mu}L(\\mu, \\rho) = 0 \\quad&\\Rightarrow\\quad \\nabla_{\\mu}\\left( -w^{T}\\mu - \\rho\\big( \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\big) \\right) = 0, \\\\\n", + " c(\\mu) \\geq 0 \\quad&\\Rightarrow\\quad \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\geq 0, \\\\\n", + " \\rho \\geq 0 \\quad&\\Rightarrow\\quad \\rho\\geq0, \\\\\n", + " \\rho c(\\mu) = 0 \\quad&\\Rightarrow\\quad \\rho\\left(\\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\right) = 0.\n", + "\\end{align}\n", + "From (1) we have that\n", + "$$\n", + " -w + \\rho 2\\Omega^{-1}(\\mu-\\bar{\\mu}) = 0 \\quad\\Rightarrow\\quad \\mu - \\bar{\\mu} = \\frac{1}{2\\rho}\\Omega w \\qquad (5)\n", + "$$\n", + "which when substituted into (4) gives\n", + "\\begin{align*}\n", + " &\\phantom{\\Rightarrow}\\quad \\rho\\left( \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\right) = 0 \\\\\n", + " &\\Rightarrow\\quad \\rho\\left( \\kappa^2 - \\big(\\frac{1}{2\\rho}\\Omega w\\big)^{T}\\Omega^{-1}\\big(\\frac{1}{2\\rho}\\Omega w\\big) \\right) = 0 \\\\\n", + " &\\Rightarrow\\quad \\rho\\kappa^2 - \\frac{1}{4\\rho}w^{T}\\Omega^{T}\\Omega^{-1}\\Omega w = 0 \\\\\n", + " &\\Rightarrow\\quad \\rho^2\\kappa^2 = \\frac{1}{4}w^{T}\\Omega w\\quad \\text{(since $\\Omega$ symmetric)} \\\\\n", + " &\\Rightarrow\\quad \\rho\\kappa = \\frac{1}{2}\\sqrt{w^{T}\\Omega w}\\quad \\text{(since $\\rho,\\kappa\\geq0$)} \\\\\n", + " &\\Rightarrow\\quad \\frac{1}{2\\rho} = \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}.\n", + "\\end{align*}\n", + "Substituting this back into (5), we find that\n", + "$$\n", + " \\mu - \\bar{\\mu} = \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w \\quad\\Rightarrow\\quad \\mu = \\bar{\\mu} + \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w\n", + "$$\n", + "as the solution for the inner problem. Substituting this back into the outer problem gives\n", + "$$\n", + " \\min_w \\left(\\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\left( \\bar{\\mu} + \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w \\right)\\right)\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + "$$\n", + "which after simplifying down and rationalizing the denominator becomes\n", + "$$\n", + " \\min_w \\left(\\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\bar{\\mu} - \\lambda\\kappa\\sqrt{w^{T}\\Omega w}\\right)\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u\n", + "$$\n", + "where $\\kappa\\in\\mathbb{R}$ is our robust optimization parameter. Since our objective has gained an additional square root term, this is obviously no longer a quadratic problem and `qpsolvers` is no longer a viable tool. We will instead now need to work with the Gurobi API directly." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Using Gurobi\n", + "\n", + "To illustrate how the setup changes when uncertainty is added, we will first look at how Gurobi handles the standard problem. The following code returns the same solution as " + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n", + "\n", + "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", + "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", + "\n", + "Optimize a model with 22 rows, 3 columns and 24 nonzeros\n", + "Model fingerprint: 0xd96d3c43\n", + "Model has 3 quadratic objective terms\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 1e+00]\n", + " Objective range [5e-01, 2e+00]\n", + " QObjective range [1e+00, 5e+00]\n", + " Bounds range [0e+00, 0e+00]\n", + " RHS range [5e-01, 1e+00]\n", + "Presolve removed 21 rows and 1 columns\n", + "Presolve time: 0.01s\n", + "Presolved: 1 rows, 2 columns, 2 nonzeros\n", + "Presolved model has 2 quadratic objective terms\n", + "Ordering time: 0.00s\n", + "\n", + "Barrier statistics:\n", + " AA' NZ : 0.000e+00\n", + " Factor NZ : 1.000e+00\n", + " Factor Ops : 1.000e+00 (less than 1 second per iteration)\n", + " Threads : 1\n", + "\n", + " Objective Residual\n", + "Iter Primal Dual Primal Dual Compl Time\n", + " 0 5.30849490e+05 -5.31975390e+05 1.50e+03 0.00e+00 3.13e+05 0s\n", + " 1 7.57990696e+02 -1.61071935e+03 3.58e+01 1.71e-13 7.91e+03 0s\n", + " 2 -7.64387653e-01 -8.13186089e+02 3.58e-05 1.23e-14 2.03e+02 0s\n", + " 3 -7.64406933e-01 -2.00154147e+00 1.88e-08 6.43e-14 3.09e-01 0s\n", + " 4 -7.99897456e-01 -9.97000994e-01 1.46e-09 4.86e-15 4.93e-02 0s\n", + " 5 -8.12408202e-01 -8.31800574e-01 1.33e-15 2.78e-17 4.85e-03 0s\n", + " 6 -8.12499997e-01 -8.12540276e-01 0.00e+00 0.00e+00 1.01e-05 0s\n", + " 7 -8.12500000e-01 -8.12500040e-01 0.00e+00 2.78e-17 1.01e-08 0s\n", + " 8 -8.12500000e-01 -8.12500000e-01 2.22e-16 5.55e-17 1.01e-11 0s\n", + "\n", + "Barrier solved model in 8 iterations and 0.02 seconds (0.00 work units)\n", + "Optimal objective -8.12500000e-01\n", + "\n", + "w = [0.50000 0.37500 0.12500]\n" + ] + } + ], + "source": [ + "# create a model for standard genetic selection\n", + "model = gp.Model(\"standardGS\")\n", + "\n", + "# define variable of interest as a continuous \n", + "w = model.addMVar(shape=problem_size, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "\n", + "# set the objective function\n", + "model.setObjective(\n", + " 0.5*w.transpose()@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values,\n", + "GRB.MINIMIZE)\n", + "\n", + "# add sub-to-half constraints\n", + "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", + "# add weight-bound constraints\n", + "model.addConstr(w >= lower_bound, name=\"lower bound\")\n", + "model.addConstr(w <= upper_bound, name=\"upper bound\")\n", + "\n", + "# solve the problem with Gurobi\n", + "model.optimize()\n", + "print(f\"w = {w.X}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Unfortunately Gurobi cannot handle our problem with its objective function in the form\n", + "$$\n", + " \\min_w \\left(\\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\lambda\\kappa\\sqrt{w^{T}\\Omega w}\\right)\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u,\n", + "$$\n", + "so further adjustments are needed first. If we define a real auxillary variable $z\\geq0$ such that $z\\leq\\sqrt{w^{T}\\Omega w}$, then our problem becomes\n", + "$$\n", + " \\min_w \\left(\\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\lambda\\kappa z\\right)\\ \\text{ s.t. }\\ z\\leq\\sqrt{w^{T}\\Omega w},\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + "$$\n", + "Since $\\kappa,\\lambda\\geq0$, and we're minimizing an objective containing \"$-\\lambda\\kappa z$\", this term of the objective will be smallest when $z>0$ is biggest. This happens precisely when it attains its upper bound from the constraint $z\\leq\\sqrt{w^{T}\\Omega w}$, hence our relaxation is justified since $z$ will push upwards.\n", + "\n", + "However, Gurobi _still_ can't handle this due to the presence of the square root, so we further make use of both $z$ and $\\sqrt{w^{T}\\Omega w}$ being positive to note that $z\\leq\\sqrt{w^{T}\\Omega w}$ can be squared on both sides:\n", + "$$\n", + " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\lambda\\kappa z\\ \\text{ s.t. }\\ z^2\\leq w^{T}\\Omega w,\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + "$$" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "To explore this with our toy problem from before, we will define\n", + "$$\n", + " \\bar{\\mu} = \\begin{bmatrix} 1 \\\\ 5 \\\\ 2 \\end{bmatrix},\\quad\n", + " \\Omega = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 4 & 0 \\\\ 0 & 0 & \\frac{1}{8} \\end{bmatrix},\\quad\n", + " \\kappa = 0.5\n", + "$$\n", + "and retain the other variables as before." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "omega = np.array([\n", + " [1, 0, 0],\n", + " [0, 4, 0],\n", + " [0, 0, 1/8]\n", + "])\n", + "\n", + "kappa = 0.5" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We then formulate this in Python as follows." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n", + "\n", + "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", + "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", + "\n", + "Optimize a model with 22 rows, 4 columns and 24 nonzeros\n", + "Model fingerprint: 0xf7819a35\n", + "Model has 3 quadratic objective terms\n", + "Model has 1 quadratic constraint\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 1e+00]\n", + " QMatrix range [1e-01, 4e+00]\n", + " Objective range [2e-01, 2e+00]\n", + " QObjective range [1e+00, 5e+00]\n", + " Bounds range [0e+00, 0e+00]\n", + " RHS range [5e-01, 1e+00]\n", + "Presolve removed 21 rows and 1 columns\n", + "\n", + "Continuous model is non-convex -- solving as a MIP\n", + "\n", + "Presolve removed 21 rows and 1 columns\n", + "Presolve time: 0.00s\n", + "Presolved: 4 rows, 6 columns, 9 nonzeros\n", + "Presolved model has 2 quadratic objective terms\n", + "Presolved model has 1 quadratic constraint(s)\n", + "Presolved model has 2 bilinear constraint(s)\n", + "Variable types: 6 continuous, 0 integer (0 binary)\n", + "Found heuristic solution: objective -1.0491801\n", + "\n", + "Root relaxation: objective -1.095481e+00, 4 iterations, 0.00 seconds (0.00 work units)\n", + "\n", + " Nodes | Current Node | Objective Bounds | Work\n", + " Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n", + "\n", + " 0 0 -1.09548 0 - -1.04918 -1.09548 4.41% - 0s\n", + " 0 0 -1.09548 0 1 -1.04918 -1.09548 4.41% - 0s\n", + " 0 0 -1.08932 0 2 -1.04918 -1.08932 3.83% - 0s\n", + " 0 0 -1.08772 0 2 -1.04918 -1.08772 3.67% - 0s\n", + " 0 0 -1.08664 0 2 -1.04918 -1.08664 3.57% - 0s\n", + " 0 0 -1.08587 0 2 -1.04918 -1.08587 3.50% - 0s\n", + " 0 0 -1.08527 0 1 -1.04918 -1.08527 3.44% - 0s\n", + " 0 0 -1.08481 0 1 -1.04918 -1.08481 3.40% - 0s\n", + " 0 0 -1.08443 0 1 -1.04918 -1.08443 3.36% - 0s\n", + " 0 0 -1.08415 0 2 -1.04918 -1.08415 3.33% - 0s\n", + " 0 0 -1.04960 0 2 -1.04918 -1.04960 0.04% - 0s\n", + " 0 0 -1.04955 0 2 -1.04918 -1.04955 0.04% - 0s\n", + " 0 2 -1.04955 0 2 -1.04918 -1.04955 0.04% - 0s\n", + "\n", + "Explored 3 nodes (22 simplex iterations) in 0.04 seconds (0.00 work units)\n", + "Thread count was 8 (of 8 available processors)\n", + "\n", + "Solution count 1: -1.04918 \n", + "\n", + "Optimal solution found (tolerance 1.00e-04)\n", + "Best objective -1.049180143116e+00, best bound -1.049222787703e+00, gap 0.0041%\n", + "w = [0.50000 0.42869 0.07131],\n", + "z = 0.9928459846964108.\n" + ] + } + ], + "source": [ + "# create a new model for robust genetic selection\n", + "model = gp.Model(\"robustGS\")\n", + "\n", + "# define variables of interest as a continuous\n", + "w = model.addMVar(shape=problem_size, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "z = model.addVar(lb=0.0, name=\"z\")\n", + "\n", + "# setup the robust objective function\n", + "model.setObjective(\n", + " 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values - lam*kappa*z,\n", + "GRB.MINIMIZE)\n", + "\n", + "# add quadratic uncertainty constraint\n", + "model.addConstr(z**2 <= w.transpose()@omega@w, name=\"uncertainty\")\n", + "# add sub-to-half constraints\n", + "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", + "# add weight-bound constraints\n", + "model.addConstr(w >= lower_bound, name=\"lower bound\")\n", + "model.addConstr(w <= upper_bound, name=\"upper bound\")\n", + "\n", + "# solve the problem with Gurobi\n", + "model.optimize()\n", + "print(f\"w = {w.X},\\nz = {z.X}.\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can repeat this experiment with varying $\\kappa$ to see how how different tolerances for uncertainty impact the robust turning point returned.\n", + "\n", + "| $\\kappa$ | $w$ | $z$ | $f(w,z)$ | Gap |\n", + "| ---: | ------------------------: | ------: | -------: | --: |\n", + "| 0.0 | [0.50000 0.37500 0.12500] | 0.0 | -8.12500000e-01 | 0% |\n", + "| 0.5 | [0.50000 0.42869 0.07131] | 0.9928459846964108 | -1.049180143116e+00 | 0.0041% |\n", + "| 1.0 | [0.50000 0.48606 0.01394] | 1.0931753188733573 | -1.309752491846e+00 | 0.0023% |\n", + "| 2.0 | [0.50000 0.50000 0.00000] | 1.1180339796124110 | -1.868033983797e+00 | 0.0005% |\n", + "| 4.0 | [0.50000 0.50000 0.00000] | 1.1180339856336740 | -2.986067972548e+00 | 0.0021%\n", + "\n", + "We can see that in the case $\\kappa = 0$ we return the standard optimization solution, as expected. We also note that there are not further changes once $\\kappa$ falls above a certain threshold, indicating that we have successfully converged to the \"riskiest\" portfolio." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Real (Simulated) Data\n", + "\n", + "Now we've looked at how these problems can be approached with Gurobi, we can try an example with realistic simulated data. Here we have data for a cohort of 50 candidates split across three data files, each containing space-separated values.\n", + "\n", + "1. In `A50.txt` the matrix $\\Sigma$ is described, with columns for $i$, $j$, and $a_{ij}$, where only the upper triangle is described since the matrix is symmetric.\n", + "2. In `EBV50.txt` the vector $\\bar{\\mu}$ is described, with only one column containing the posterior mean over 1000 samples.\n", + "3. In `S50.txt` the matrix $\\Omega$ is described, with columns for $i$, $j$, and $\\sigma_{ij}$, where only the upper triangle is described since the matrix is symmetric.\n", + "\n", + "Note that this particular problem does not contain a separate file for sex data; odd indexed candidates are sires and even indexed candidates are dams. For now, we also don't have a $l$ or $u$ bounding the possible weights.\n", + "\n", + "\n", + "The first step then is clearly going to be loading the matrices in particular from file. We create the following function to do so.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "def load_problem(A_filename, E_filename, S_filename, dimension=False):\n", + " \"\"\"\n", + " Used to load genetic selection problems into NumPy. It takes three\n", + " string inputs for filenames where Sigma, Mu, and Omega are stored,\n", + " as well as an optional integer input for problem dimension if this\n", + " is known. If it's know know, it's worked out based on E_filename.\n", + "\n", + " As output, it returns (A, E, S, n), where A and S are n-by-n NumPy\n", + " arrays, E is a length n NumPy array, and n is an integer.\n", + " \"\"\"\n", + "\n", + " def load_symmetric_matrix(filename, dimension):\n", + " \"\"\"\n", + " Since NumPy doesn't have a stock way to load matrices\n", + " stored in coordinate format format, this adds one.\n", + " \"\"\"\n", + "\n", + " matrix = np.zeros([dimension, dimension])\n", + "\n", + " with open(filename, 'r') as file:\n", + " for line in file:\n", + " i, j, entry = line.split(\" \")\n", + " # data files indexed from 1, not 0\n", + " matrix[int(i)-1, int(j)-1] = entry\n", + " matrix[int(j)-1, int(i)-1] = entry\n", + "\n", + " return matrix\n", + "\n", + "\n", + " # if dimension wasn't supplied, need to find that\n", + " if not dimension:\n", + " # get dimension from EBV, since it's the smallest file\n", + " with open(E_filename, 'r') as file:\n", + " dimension = sum(1 for _ in file)\n", + "\n", + " # EBV isn't in coordinate format so can be loaded directly\n", + " E = np.loadtxt(E_filename) \n", + " # A and S are stored by coordinates so need special loader\n", + " A = load_symmetric_matrix(A_filename, dimension)\n", + " S = load_symmetric_matrix(S_filename, dimension)\n", + "\n", + " return A, E, S, dimension" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For the given example problem we can now solve it with Gurobi using the exact same methods as before, for both the standard genetic selection problem and the robust version of the problem. The following cell does both alongside each other to accentuate the differences." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Set parameter TimeLimit to value 300\n", + "Set parameter MIPGap to value 0.05\n", + "Set parameter MIPGap to value 0.05\n", + "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n", + "\n", + "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", + "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", + "\n", + "Optimize a model with 4 rows, 50 columns and 100 nonzeros\n", + "Model fingerprint: 0x534ff85e\n", + "Model has 1275 quadratic objective terms\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 1e+00]\n", + " Objective range [7e-01, 1e+00]\n", + " QObjective range [5e-02, 1e+00]\n", + " Bounds range [0e+00, 0e+00]\n", + " RHS range [5e-01, 5e-01]\n", + "Presolve removed 2 rows and 0 columns\n", + "Presolve time: 0.01s\n", + "Presolved: 2 rows, 50 columns, 50 nonzeros\n", + "Presolved model has 1275 quadratic objective terms\n", + "Ordering time: 0.00s\n", + "\n", + "Barrier statistics:\n", + " Free vars : 49\n", + " AA' NZ : 1.274e+03\n", + " Factor NZ : 1.326e+03\n", + " Factor Ops : 4.553e+04 (less than 1 second per iteration)\n", + " Threads : 1\n", + "\n", + " Objective Residual\n", + "Iter Primal Dual Primal Dual Compl Time\n", + " 0 3.76397773e+05 -4.22619805e+05 2.50e+04 1.12e+00 9.99e+05 0s\n", + " 1 -5.01117677e+01 -9.99230833e+02 2.68e+01 1.25e-04 1.09e+03 0s\n", + " 2 -8.81632395e-01 -9.63667643e+02 2.68e-05 1.25e-10 1.93e+01 0s\n", + " 3 -8.81606636e-01 -2.03034027e+00 5.19e-09 2.41e-14 2.30e-02 0s\n", + " 4 -9.04332054e-01 -1.04180295e+00 4.73e-10 2.19e-15 2.75e-03 0s\n", + " 5 -9.59656693e-01 -1.01318293e+00 1.67e-15 8.88e-16 1.07e-03 0s\n", + " 6 -9.73793052e-01 -9.76849441e-01 2.22e-16 1.33e-15 6.11e-05 0s\n", + " 7 -9.76163966e-01 -9.76239391e-01 1.39e-15 6.66e-16 1.51e-06 0s\n", + " 8 -9.76213879e-01 -9.76217890e-01 9.91e-15 5.55e-16 8.02e-08 0s\n", + " 9 -9.76215229e-01 -9.76215592e-01 1.65e-13 5.55e-16 7.25e-09 0s\n", + " 10 -9.76215438e-01 -9.76215482e-01 7.03e-13 6.66e-16 8.86e-10 0s\n", + " 11 -9.76215473e-01 -9.76215476e-01 4.24e-13 6.66e-16 4.91e-11 0s\n", + "\n", + "Barrier solved model in 11 iterations and 0.02 seconds (0.00 work units)\n", + "Optimal objective -9.76215473e-01\n", + "\n", + "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n", + "\n", + "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", + "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", + "\n", + "Optimize a model with 4 rows, 51 columns and 100 nonzeros\n", + "Model fingerprint: 0xd47c065e\n", + "Model has 1275 quadratic objective terms\n", + "Model has 1 quadratic constraint\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 1e+00]\n", + " QMatrix range [1e-05, 1e+00]\n", + " Objective range [7e-01, 1e+00]\n", + " QObjective range [5e-02, 1e+00]\n", + " Bounds range [0e+00, 0e+00]\n", + " RHS range [5e-01, 5e-01]\n", + "Presolve removed 2 rows and 0 columns\n", + "\n", + "Continuous model is non-convex -- solving as a MIP\n", + "\n", + "Presolve removed 2 rows and 0 columns\n", + "Presolve time: 0.00s\n", + "Presolved: 2503 rows, 1327 columns, 6326 nonzeros\n", + "Presolved model has 1275 quadratic objective terms\n", + "Presolved model has 1 quadratic constraint(s)\n", + "Presolved model has 1275 bilinear constraint(s)\n", + "Variable types: 1327 continuous, 0 integer (0 binary)\n", + "Found heuristic solution: objective -1.1617295\n", + "\n", + "Root relaxation: objective -3.490954e+00, 43 iterations, 0.00 seconds (0.00 work units)\n", + "\n", + " Nodes | Current Node | Objective Bounds | Work\n", + " Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n", + "\n", + " 0 0 -3.49095 0 38 -1.16173 -3.49095 200% - 0s\n", + " 0 0 -3.49095 0 42 -1.16173 -3.49095 200% - 0s\n", + " 0 0 -1.58157 0 231 -1.16173 -1.58157 36.1% - 0s\n", + " 0 0 -1.34540 0 512 -1.16173 -1.34540 15.8% - 0s\n", + " 0 0 -1.32278 0 619 -1.16173 -1.32278 13.9% - 0s\n", + " 0 0 -1.32278 0 59 -1.16173 -1.32278 13.9% - 1s\n", + " 0 0 -1.32278 0 59 -1.16173 -1.32278 13.9% - 1s\n", + " 0 2 -1.32278 0 59 -1.16173 -1.32278 13.9% - 1s\n", + " 110 118 -1.23808 14 191 -1.16173 -1.27940 10.1% 343 5s\n", + " 340 353 -1.19438 40 319 -1.16173 -1.27940 10.1% 377 10s\n", + "H 551 573 -1.1617295 -1.27240 9.53% 438 12s\n", + "\n", + "Explored 600 nodes (256270 simplex iterations) in 13.46 seconds (12.04 work units)\n", + "Thread count was 8 (of 8 available processors)\n", + "\n", + "Solution count 1: -1.16173 \n", + "\n", + "Optimal solution found (tolerance 5.00e-02)\n", + "Best objective -1.161729466120e+00, best bound -1.217554010902e+00, gap 4.8053%\n", + "\n", + "SIRE WEIGHTS\t\t\t DAM WEIGHTS\n", + "--------------------\t\t --------------------\n", + " i w_std w_rbs\t\t i w_std w_rbs\n", + "00 0.00000 0.00000 01 0.08319 0.14976\n", + "02 0.11649 0.49337 03 0.00000 0.00000\n", + "04 0.00000 0.00000 05 0.00000 0.00000\n", + "06 0.00000 0.00000 07 0.06817 0.06251\n", + "08 0.00000 0.00000 09 0.03762 0.00000\n", + "10 0.00000 0.00000 11 0.00000 0.00000\n", + "12 0.00000 0.00000 13 0.00355 0.00000\n", + "14 0.11542 0.00518 15 0.06966 0.00000\n", + "16 0.00000 0.00000 17 0.00000 0.00000\n", + "18 0.00000 0.00000 19 0.00000 0.00000\n", + "20 0.00000 0.00000 21 0.00001 0.00000\n", + "22 0.00000 0.00000 23 0.00000 0.00000\n", + "24 0.03537 0.00000 25 0.00000 0.00000\n", + "26 0.00000 0.00000 27 0.00000 0.00000\n", + "28 0.00000 0.00000 29 0.00000 0.00000\n", + "30 0.00000 0.00000 31 0.03580 0.00000\n", + "32 0.08447 0.00000 33 0.00000 0.00000\n", + "34 0.00000 0.00000 35 0.02849 0.09487\n", + "36 0.13489 0.00144 37 0.07446 0.10716\n", + "38 0.00000 0.00000 39 0.00000 0.00000\n", + "40 0.00000 0.00000 41 0.07312 0.08570\n", + "42 0.00000 0.00000 43 0.00000 0.00000\n", + "44 0.01336 0.00000 45 0.00000 0.00000\n", + "46 0.00000 0.00000 47 0.02593 0.00000\n", + "48 0.00000 0.00000 49 0.00000 0.00000\n", + "\n", + "Maximum change: 0.37689\n", + "Average change: 0.02220\n", + "Minimum change: 0.00000\n" + ] + } + ], + "source": [ + "sigma, mubar, omega, n = load_problem(\n", + " \"../Example/A50.txt\",\n", + " \"../Example/EBV50.txt\",\n", + " \"../Example/S50.txt\",\n", + " 50)\n", + "\n", + "lam = 0.5\n", + "kappa = 2\n", + "\n", + "# define the M so that column i is [1;0] if i is a sire (so even) and [0;1] otherwise \n", + "M = np.zeros((2, n))\n", + "M[0, range(0,50,2)] = 1\n", + "M[1, range(1,50,2)] = 1\n", + "# define the right hand side of the constraint Mx = m\n", + "m = np.array([[0.5], [0.5]])\n", + "\n", + "# create models for standard and robust genetic selection\n", + "model_std = gp.Model(\"n50standard\")\n", + "model_rbs = gp.Model(\"n50robust\")\n", + "\n", + "# initialise w for both models, z for robust model\n", + "w_std = model_std.addMVar(shape=n, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\") \n", + "w_rbs = model_rbs.addMVar(shape=n, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "z_rbs = model_rbs.addVar(lb=0.0, name=\"z\")\n", + "\n", + "# define the objective functions for both models\n", + "model_std.setObjective(\n", + " 0.5*w_std.transpose()@(sigma@w_std) - lam*w_std.transpose()@mubar,\n", + "GRB.MINIMIZE)\n", + "\n", + "model_rbs.setObjective(\n", + " # Gurobi does offer a way to set one objective in terms of another, i.e.\n", + " # we could use `model_std.getObjective() - lam*kappa*z_rbs` to define this\n", + " # robust objective, but it results in a significant slowdown in code.\n", + " 0.5*w_rbs.transpose()@(sigma@w_rbs) - lam*w_rbs.transpose()@mubar - lam*kappa*z_rbs,\n", + "GRB.MINIMIZE)\n", + "\n", + "# add sum-to-half constraints to both models\n", + "model_std.addConstr(M @ w_std == m, name=\"sum-to-half\")\n", + "model_rbs.addConstr(M @ w_rbs == m, name=\"sum-to-half\")\n", + "\n", + "# add quadratic uncertainty constraint to the robust model\n", + "model_rbs.addConstr(z_rbs**2 <= w_rbs.transpose()@omega@w_rbs, name=\"uncertainty\")\n", + "\n", + "# since working with non-trivial size, set a time limit\n", + "time_limit = 60*5 # 5 minutes\n", + "model_std.setParam(GRB.Param.TimeLimit, time_limit)\n", + "model_std.setParam(GRB.Param.TimeLimit, time_limit)\n", + "\n", + "# for the same reason, also set a duality gap tolerance\n", + "duality_gap = 0.05\n", + "model_std.setParam('MIPGap', duality_gap)\n", + "model_rbs.setParam('MIPGap', duality_gap)\n", + "\n", + "# solve both problems with Gurobi\n", + "model_std.optimize()\n", + "model_rbs.optimize()\n", + "\n", + "# HACK code which prints the results for comparison in a nice format\n", + "print(\"\\nSIRE WEIGHTS\\t\\t\\t DAM WEIGHTS\")\n", + "print(\"-\"*20 + \"\\t\\t \" + \"-\"*20)\n", + "print(\" i w_std w_rbs\\t\\t i w_std w_rbs\")\n", + "for candidate in range(25):\n", + " print(f\"{candidate*2:02d} {w_std.X[candidate*2]:.5f} {w_rbs.X[candidate*2]:.5f} \\\n", + " {candidate*2+1:02d} {w_std.X[candidate*2+1]:.5f} {w_rbs.X[candidate*2+1]:.5f}\")\n", + " \n", + "print(f\"\\nMaximum change: {max(np.abs(w_std.X-w_rbs.X)):.5f}\")\n", + "print(f\"Average change: {np.mean(np.abs(w_std.X-w_rbs.X)):.5f}\")\n", + "print(f\"Minimum change: {min(np.abs(w_std.X-w_rbs.X)):.5f}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can see that with parameters $\\lambda = 0.5, \\kappa = 2$, there's a more significant shift in the portfolios produced and we see in the return maximization version of the problem. It also takes Gurobi significantly longer to compute. More testing would be needed to confirm why this was the case." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/robust-genetics.ipynb b/robust-genetics.ipynb index 178193f..1507c41 100644 --- a/robust-genetics.ipynb +++ b/robust-genetics.ipynb @@ -57,13 +57,25 @@ "source": [ "## Standard Problem\n", "\n", - "In the context of genetic selection, we want to maximize selection of genetic merit (a measure of desirable traits) while minimizing risks due to inbreeding. This can be formed mathematically as\n", + "In the context of genetic selection, we want to maximize selection of genetic merit (a measure of desirable traits) while minimizing risks due to inbreeding. This can be formed mathematically as a multi-objective optimization problem\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu\\ \\text{ subject to }\\ w_{\\mathcal{S}}^{T}e_{\\mathcal{S}}^{} = \\frac{1}{2},\\ w_{\\mathcal{D}}^{T}e_{\\mathcal{D}}^{} = \\frac{1}{2},\\ l\\leq w\\leq u,\n", + " \\min_w \\left(\\frac{1}{2}w^{T}\\Sigma w\\right),\\ \\max_w \\left(w^{T}\\mu\\right) \\text{ subject to }\\ w_{\\mathcal{S}}^{T}e_{\\mathcal{S}}^{} = \\frac{1}{2},\\ w_{\\mathcal{D}}^{T}e_{\\mathcal{D}}^{} = \\frac{1}{2},\\ l\\leq w\\leq u\n", "$$\n", "where $w$ is the vector of proportional contributions, $\\Sigma$ is a matrix encoding risk, $\\mu$ is a vector encoding returns, $l$ encodes lower bounds on contributions, $u$ encodes upper bounds on contributions, $\\mathcal{S}$ is an index set of candidates who are sires, and $\\mathcal{D}$ is an index set of candidates who are dams.\n", "\n", - "In this representation of the problem, $\\lambda$ is a control variable which balances how we trade of between risk and return. Each value of $\\lambda$ will give a different solution on the critical frontier of the problem." + "When using Critical Line Algorithm we typically frame the problem as minimizing risk for a target level of return, $\\mu_p$, and can formulate the problem as a single objective\n", + "$$\n", + " \\min_w \\left(\\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu\\right) \\text{ subject to }\\ w_{\\mathcal{S}}^{T}e_{\\mathcal{S}}^{} = \\frac{1}{2},\\ w_{\\mathcal{D}}^{T}e_{\\mathcal{D}}^{} = \\frac{1}{2},\\ l\\leq w\\leq u.\n", + "$$\n", + "In this representation of the problem, $\\lambda$ is a control variable which balances how we trade of between risk and return. Each value of $\\lambda$ will give a different solution on the critical frontier of the problem.\n", + "\n", + "The mathematics of framing it this way works out nicer than maximizing return for a target level of risk, $\\sigma_p$. However, [we will find](Scratch/robust-genetics-reframe.ipynb) that this framing actually makes the mathematics _more_ complicated in the context of robust optimization.\n", + "\n", + "For that reason we will instead reframe the problem as \n", + "$$\n", + " \\max_w \\left(w^{T}\\mu - \\frac{\\lambda}{2}w^{T}\\Sigma w\\right) \\text{ subject to }\\ w_{\\mathcal{S}}^{T}e_{\\mathcal{S}}^{} = \\frac{1}{2},\\ w_{\\mathcal{D}}^{T}e_{\\mathcal{D}}^{} = \\frac{1}{2},\\ l\\leq w\\leq u.\n", + "$$\n", + "Note that the $\\lambda$ here corresponding to a particular point on the frontier will not have the same value as the $\\lambda$ for the same point in the other framing." ] }, { @@ -97,7 +109,7 @@ "source": [ "### Toy Example ($n = 3$)\n", "\n", - "Lets see how this works in an example. We will start by looking how this problem might be solving using Python's [qpsolvers](https://qpsolvers.github.io/qpsolvers/index.html) library. Consider the problem where\n", + "Lets see how this works in an example. We will start by looking how this problem might be solving using a higher level modelling package, Python's [qpsolvers](https://qpsolvers.github.io/qpsolvers/index.html) library. Consider the problem where\n", "$$\n", " \\mu = \\begin{bmatrix} 1 \\\\ 5 \\\\ 2 \\end{bmatrix},\\quad\n", " \\Sigma = \\begin{bmatrix} 1 & 0 & 0 \\\\ 0 & 5 & 0 \\\\ 0 & 0 & 3 \\end{bmatrix}, \\quad\n", @@ -147,7 +159,7 @@ "outputs": [], "source": [ "# OPTIMIZATION SETUP VARIABLES\n", - "lam = 0.5\n", + "lam = 2\n", "# define the M so that column i is [1;0] if i is a sire and [0;1] otherwise \n", "M = np.zeros((2, problem_size))\n", "M[0, sire_indices] = 1\n", @@ -182,8 +194,8 @@ "# SOLVE THE PROBLEM\n", "def solveGS(lam):\n", " return solve_qp(\n", - " P = relationship_matrix,\n", - " q = -lam*expected_breeding_values,\n", + " P = lam*relationship_matrix,\n", + " q = -expected_breeding_values,\n", " G = None,\n", " h = None,\n", " A = M,\n", @@ -200,7 +212,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Excellent, we have code which is telling us that our optimal contributions for $\\lambda = 0.5$ are $w_1 = 0.5$, $w_2 = 0.375$, and $w_3 = 0.125$. We could vary our $\\lambda$ value too and find other points on the frontier." + "Excellent, we have code which is telling us that our optimal contributions for $\\lambda = 2$ are $w_1 = 0.5$, $w_2 = 0.375$, and $w_3 = 0.125$. We could vary our $\\lambda$ value too and find other points on the frontier.¹" ] }, { @@ -212,22 +224,17 @@ "name": "stdout", "output_type": "stream", "text": [ - "lambda = 0.0: w = [0.50000 0.18750 0.31250]\n", - "lambda = 0.2: w = [0.50000 0.26250 0.23750]\n", - "lambda = 0.4: w = [0.50000 0.33750 0.16250]\n", - "lambda = 0.6: w = [0.50000 0.41250 0.08750]\n", - "lambda = 0.8: w = [0.50000 0.48750 0.01250]\n", - "lambda = 1.0: w = [0.50000 0.50000 0.00000]\n" + "lambda = 1: w = [0.50000 0.50000 0.00000]\n", + "lambda = 2: w = [0.50000 0.37500 0.12500]\n", + "lambda = 3: w = [0.50000 0.31250 0.18750]\n", + "lambda = 4: w = [0.50000 0.28125 0.21875]\n", + "lambda = 5: w = [0.50000 0.26250 0.23750]\n" ] } ], "source": [ - "print(f\"lambda = 0.0: w = {solveGS(0.0)}\")\n", - "print(f\"lambda = 0.2: w = {solveGS(0.2)}\")\n", - "print(f\"lambda = 0.4: w = {solveGS(0.4)}\")\n", - "print(f\"lambda = 0.6: w = {solveGS(0.6)}\")\n", - "print(f\"lambda = 0.8: w = {solveGS(0.8)}\")\n", - "print(f\"lambda = 1.0: w = {solveGS(1.0)}\")" + "for x in range(1,6):\n", + " print(f\"lambda = {x}: w = {solveGS(x)}\")" ] }, { @@ -236,11 +243,11 @@ "source": [ "## Robust Optimization\n", "\n", - "A wrinkle to this is that we don't typically our problem variables with certainty. There _are_ ways to define $\\Sigma$ based on relationships within the cohort (which are known and discussed more later) but $\\mu$ is an estimated variable. In particular we say $\\mu$ has a univariate normal distribution, $\\mu\\sim N(\\bar{\\mu}, \\Omega)$. This means we must turn to optimization tools which can address this uncertainty.¹\n", + "A wrinkle to this is that we don't typically our problem variables with certainty. There _are_ ways to define $\\Sigma$ based on relationships within the cohort (which are known and discussed more later) but $\\mu$ is an estimated variable. In particular we say $\\mu$ has a univariate normal distribution, $\\mu\\sim N(\\bar{\\mu}, \\Omega)$. This means we must turn to optimization tools which can address this uncertainty.²\n", "\n", "Robust optimization is one such tool in which we adjust the objective function to model the inherent uncertainty in the problem. Within this framework, the problem would be the bilevel optimization problem\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda\\left( \\max_{\\mu\\in U_{\\mu}} w^{T}\\mu \\right)\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + " \\max_w \\left(\\min_{\\mu\\in U_{\\mu}} \\left(w^{T}\\mu \\right) - \\frac{\\lambda}{2}w^{T}\\Sigma w\\right) \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", "$$\n", "The result of the inner problem determines the outer problem in a follower-leader setup. However, the inner problem is relatively simple and can be solved explicitly based on our choice of $U_{\\mu}$, representing an uncertainty set around the value of $\\mu$.\n" ] @@ -251,34 +258,34 @@ "source": [ "### Quadratic Uncertainty Sets\n", "\n", - "There are many ways to define the uncertainty set $U_{\\mu}$, but for practical reasons relating to continuity and differentiability² it's beneficial to define it using a quadratic uncertainty set. We say that\n", + "There are many ways to define the uncertainty set $U_{\\mu}$, but for practical reasons relating to continuity and differentiability it's beneficial to define it using a quadratic uncertainty set.³ We say that\n", "$$\n", " U_{\\mu} := \\left\\lbrace \\mu :\\ {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\leq \\kappa^2 \\right\\rbrace\n", "$$\n", - "which in effect bounds the uncertainty in a ball about $\\bar{\\mu}$ with parameter $\\kappa$.³ This means that our lower level problem $\\max_{\\mu\\in U_{\\mu}} w^{T}\\mu$ becomes\n", + "which in effect bounds the uncertainty in a ball about $\\bar{\\mu}$ with parameter $\\kappa$.⁴ This means that our lower level problem $\\min_{\\mu\\in U_{\\mu}} (w^{T}\\mu)$ becomes\n", "$$\n", - " \\max_{\\mu} w^{T}\\mu\\quad \\text{subject to}\\quad {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\leq \\kappa^2,\n", + " \\min_{\\mu} \\left(w^{T}\\mu \\right)\\quad \\text{subject to}\\quad {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\leq \\kappa^2,\n", "$$\n", "or, in standard form,\n", "$$\n", - " \\min_{\\mu} \\left(-w^{T}\\mu\\right) \\quad\\text{subject to}\\quad \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\geq 0.\n", + " \\min_{\\mu} \\left(w^{T}\\mu\\right) \\quad\\text{subject to}\\quad \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\geq 0.\n", "$$\n", - "This is convex, since $\\Omega$ is a positive definite matrix and $-w^{T}\\mu$ is an affine function. If we define a Lagrangian multiplier $\\rho\\in\\mathbb{R}$, the KKT conditions of this problem are:\n", + "This is convex, since $\\Omega$ is a positive definite matrix and $w^{T}\\mu$ is an affine function. If we define a Lagrangian multiplier $\\rho\\in\\mathbb{R}$, the KKT conditions of this problem are:\n", "\\begin{align}\n", - " \\nabla_{\\mu}L(\\mu, \\rho) = 0 \\quad&\\Rightarrow\\quad \\nabla_{\\mu}\\left( -w^{T}\\mu - \\rho\\big( \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\big) \\right) = 0, \\\\\n", + " \\nabla_{\\mu}L(\\mu, \\rho) = 0 \\quad&\\Rightarrow\\quad \\nabla_{\\mu}\\left( w^{T}\\mu - \\rho\\big( \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\big) \\right) = 0, \\\\\n", " c(\\mu) \\geq 0 \\quad&\\Rightarrow\\quad \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\geq 0, \\\\\n", " \\rho \\geq 0 \\quad&\\Rightarrow\\quad \\rho\\geq0, \\\\\n", " \\rho c(\\mu) = 0 \\quad&\\Rightarrow\\quad \\rho\\left(\\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\right) = 0.\n", "\\end{align}\n", "From (1) we have that\n", "$$\n", - " -w + \\rho 2\\Omega^{-1}(\\mu-\\bar{\\mu}) = 0 \\quad\\Rightarrow\\quad \\mu - \\bar{\\mu} = \\frac{1}{2\\rho}\\Omega w \\qquad (5)\n", + " w + \\rho 2\\Omega^{-1}(\\mu-\\bar{\\mu}) = 0 \\quad\\Rightarrow\\quad \\mu - \\bar{\\mu} = -\\frac{1}{2\\rho}\\Omega w \\qquad (5)\n", "$$\n", "\\end{align*}\n", "which when substituted into (4) gives\n", "\\begin{align*}\n", " &\\phantom{\\Rightarrow}\\quad \\rho\\left( \\kappa^2 - {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) \\right) = 0 \\\\\n", - " &\\Rightarrow\\quad \\rho\\left( \\kappa^2 - \\big(\\frac{1}{2\\rho}\\Omega w\\big)^{T}\\Omega^{-1}\\big(\\frac{1}{2\\rho}\\Omega w\\big) \\right) = 0 \\\\\n", + " &\\Rightarrow\\quad \\rho\\left( \\kappa^2 - \\big(-\\frac{1}{2\\rho}\\Omega w\\big)^{T}\\Omega^{-1}\\big(-\\frac{1}{2\\rho}\\Omega w\\big) \\right) = 0 \\\\\n", " &\\Rightarrow\\quad \\rho\\kappa^2 - \\frac{1}{4\\rho}w^{T}\\Omega^{T}\\Omega^{-1}\\Omega w = 0 \\\\\n", " &\\Rightarrow\\quad \\rho^2\\kappa^2 = \\frac{1}{4}w^{T}\\Omega w\\quad \\text{(since $\\Omega$ symmetric)} \\\\\n", " &\\Rightarrow\\quad \\rho\\kappa = \\frac{1}{2}\\sqrt{w^{T}\\Omega w}\\quad \\text{(since $\\rho,\\kappa\\geq0$)} \\\\\n", @@ -286,15 +293,15 @@ "\\end{align*}\n", "Substituting this back into (5), we find that\n", "$$\n", - " \\mu - \\bar{\\mu} = \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w \\quad\\Rightarrow\\quad \\mu = \\bar{\\mu} + \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w\n", + " \\mu - \\bar{\\mu} = -\\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w \\quad\\Rightarrow\\quad \\mu = \\bar{\\mu} - \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w\n", "$$\n", "as the solution for the inner problem. Substituting this back into the outer problem gives\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\left( \\bar{\\mu} + \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w \\right)\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + " \\max_w \\left(w^{T}\\left( \\bar{\\mu} - \\frac{\\kappa}{\\sqrt{w^{T}\\Omega w}}\\Omega w \\right) - \\frac{\\lambda}{2}w^{T}\\Sigma w\\right) \\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", "$$\n", "which after simplifying down and rationalizing the denominator becomes\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\bar{\\mu} - \\lambda\\kappa\\sqrt{w^{T}\\Omega w} \\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u\n", + " \\max_w \\left(w^{T}\\bar{\\mu} - \\kappa\\sqrt{w^{T}\\Omega w} - \\frac{\\lambda}{2}w^{T}\\Sigma w\\right) \\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", "$$\n", "where $\\kappa\\in\\mathbb{R}$ is our robust optimization parameter. Since our objective has gained an additional square root term, this is obviously no longer a quadratic problem and `qpsolvers` is no longer a viable tool. We will instead now need to work with the Gurobi API directly." ] @@ -323,16 +330,16 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 22 rows, 3 columns and 24 nonzeros\n", - "Model fingerprint: 0xd96d3c43\n", + "Model fingerprint: 0xd5877ae8\n", "Model has 3 quadratic objective terms\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", - " Objective range [5e-01, 2e+00]\n", - " QObjective range [1e+00, 5e+00]\n", + " Objective range [1e+00, 5e+00]\n", + " QObjective range [2e+00, 1e+01]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [5e-01, 1e+00]\n", "Presolve removed 21 rows and 1 columns\n", - "Presolve time: 0.01s\n", + "Presolve time: 0.00s\n", "Presolved: 1 rows, 2 columns, 2 nonzeros\n", "Presolved model has 2 quadratic objective terms\n", "Ordering time: 0.00s\n", @@ -345,18 +352,18 @@ "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 5.30849490e+05 -5.31975390e+05 1.50e+03 0.00e+00 3.13e+05 0s\n", - " 1 7.57990696e+02 -1.61071935e+03 3.58e+01 1.71e-13 7.91e+03 0s\n", - " 2 -7.64387653e-01 -8.13186089e+02 3.58e-05 1.23e-14 2.03e+02 0s\n", - " 3 -7.64406933e-01 -2.00154147e+00 1.88e-08 6.43e-14 3.09e-01 0s\n", - " 4 -7.99897456e-01 -9.97000994e-01 1.46e-09 4.86e-15 4.93e-02 0s\n", - " 5 -8.12408202e-01 -8.31800574e-01 1.33e-15 2.78e-17 4.85e-03 0s\n", - " 6 -8.12499997e-01 -8.12540276e-01 0.00e+00 0.00e+00 1.01e-05 0s\n", - " 7 -8.12500000e-01 -8.12500040e-01 0.00e+00 2.78e-17 1.01e-08 0s\n", - " 8 -8.12500000e-01 -8.12500000e-01 2.22e-16 5.55e-17 1.01e-11 0s\n", + " 0 -4.24839898e+06 4.25290078e+06 1.50e+03 0.00e+00 3.13e+05 0s\n", + " 1 -6.22305125e+03 9.63367781e+03 3.63e+01 0.00e+00 7.90e+03 0s\n", + " 2 1.53372163e+00 3.34533125e+03 3.63e-05 0.00e+00 1.05e+02 0s\n", + " 3 1.53373576e+00 5.70008844e+00 8.94e-09 0.00e+00 1.30e-01 0s\n", + " 4 1.57092842e+00 2.18994634e+00 1.18e-09 0.00e+00 1.93e-02 0s\n", + " 5 1.61185552e+00 1.85809165e+00 1.22e-15 0.00e+00 7.69e-03 0s\n", + " 6 1.62487212e+00 1.63327073e+00 0.00e+00 0.00e+00 2.62e-04 0s\n", + " 7 1.62500000e+00 1.62501541e+00 0.00e+00 5.55e-17 4.81e-07 0s\n", + " 8 1.62500000e+00 1.62500002e+00 0.00e+00 0.00e+00 4.81e-10 0s\n", "\n", - "Barrier solved model in 8 iterations and 0.03 seconds (0.00 work units)\n", - "Optimal objective -8.12500000e-01\n", + "Barrier solved model in 8 iterations and 0.01 seconds (0.00 work units)\n", + "Optimal objective 1.62500000e+00\n", "\n", "w = [0.50000 0.37500 0.12500]\n" ] @@ -367,18 +374,18 @@ "model = gp.Model(\"standardGS\")\n", "\n", "# define variable of interest as a continuous \n", - "w = model.addMVar(shape=problem_size, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "w = model.addMVar(shape=problem_size, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\")\n", "\n", "# set the objective function\n", "model.setObjective(\n", - " 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values,\n", - "GRB.MINIMIZE)\n", + " w.transpose()@expected_breeding_values - (lam/2)*w@(relationship_matrix@w),\n", + "GRB.MAXIMIZE)\n", "\n", "# add sub-to-half constraints\n", "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", "# add weight-bound constraints\n", - "model.addConstr(w >= lower_bound, name=\"lower bound\")\n", - "model.addConstr(w <= upper_bound, name=\"upper bound\")\n", + "model.addConstr(w >= lower_bound, name=\"lower-bound\")\n", + "model.addConstr(w <= upper_bound, name=\"upper-bound\")\n", "\n", "# solve the problem with Gurobi\n", "model.optimize()\n", @@ -389,19 +396,23 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Unfortunately Gurobi cannot handle our problem with its objective function in the form\n", + "Unfortunately Gurobi cannot handle our robust problem with its objective function in the form\n", + "$$\n", + " \\max_w \\left(w^{T}\\bar{\\mu} - \\kappa\\sqrt{w^{T}\\Omega w} - \\frac{\\lambda}{2}w^{T}\\Sigma w\\right) \\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u,\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\lambda\\kappa\\sqrt{w^{T}\\Omega w}\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u,\n", + "or in its norm form\n", "$$\n", - "so some further adjustments are needed first. If we define a real auxillary variable $z\\geq0$ such that $z\\leq\\sqrt{w^{T}\\Omega w}$, then our problem becomes\n", + " \\max_w \\left(w^{T}\\bar{\\mu} - \\kappa\\|\\Omega^{\\frac{1}{2}}w\\| - \\frac{\\lambda}{2}w^{T}\\Sigma w\\right) \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u,\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\lambda\\kappa z\\ \\text{ s.t. }\\ z\\leq\\sqrt{w^{T}\\Omega w},\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + "so some further adjustments are needed first. If we define a real auxillary variable $z\\geq0$ such that $z\\geq\\sqrt{w^{T}\\Omega w}$, then our problem becomes\n", "$$\n", - "TODO: _better explain the idea of $z$ pushing up so that this switch doesn't have any practical difference._\n", + " \\max_w \\left(w^{T}\\bar{\\mu} - \\kappa z - \\frac{\\lambda}{2}w^{T}\\Sigma w\\right) \\ \\text{ subject to }\\ z\\geq\\sqrt{w^{T}\\Omega w},\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + "$$\n", + "Since $z>0$, $\\kappa\\geq0$, and we're maximizing an objective containing \"$-\\kappa z$\", this term of the objective will be biggest when $z$ is smallest. This will happen precisely when it attains its lower bound from the constraint $z\\geq\\sqrt{w^{T}\\Omega w}$, hence our relaxation is justified since $z$ will push downwards.\n", "\n", - "However, Gurobi _still_ can't handle this due to the presence of the square root, so we further make use of both $z$ and $\\sqrt{w^{T}\\Omega w}$ being positive to note that $z\\leq\\sqrt{w^{T}\\Omega w}$ can be squared on both sides:\n", + "However, Gurobi _still_ can't handle this due to the presence of the square root, so we further make use of both $z$ and $\\sqrt{w^{T}\\Omega w}$ being positive to note that $z\\geq\\sqrt{w^{T}\\Omega w}$ can be squared on both sides:\n", "$$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\lambda\\kappa z\\ \\text{ s.t. }\\ z^2\\leq w^{T}\\Omega w,\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + " \\max_w \\left(w^{T}\\bar{\\mu} - \\kappa z - \\frac{\\lambda}{2}w^{T}\\Sigma w\\right) \\ \\text{ subject to }\\ z^2\\geq w^{T}\\Omega w,\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", "$$" ] }, @@ -454,44 +465,43 @@ "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", - "Optimize a model with 23 rows, 4 columns and 25 nonzeros\n", - "Model fingerprint: 0xb99edb8a\n", + "Optimize a model with 22 rows, 4 columns and 24 nonzeros\n", + "Model fingerprint: 0xe9e68b83\n", "Model has 3 quadratic objective terms\n", - "Model has 3 quadratic constraints\n", + "Model has 1 quadratic constraint\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " QMatrix range [1e-01, 4e+00]\n", - " Objective range [2e-01, 2e+00]\n", - " QObjective range [1e+00, 5e+00]\n", + " Objective range [5e-01, 5e+00]\n", + " QObjective range [2e+00, 1e+01]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [5e-01, 1e+00]\n", - "Presolve removed 22 rows and 1 columns\n", + "Presolve removed 21 rows and 1 columns\n", "Presolve time: 0.01s\n", - "Presolved: 9 rows, 8 columns, 14 nonzeros\n", - "Presolved model has 3 second-order cone constraints\n", + "Presolved: 7 rows, 8 columns, 11 nonzeros\n", + "Presolved model has 2 second-order cone constraints\n", "Ordering time: 0.00s\n", "\n", "Barrier statistics:\n", - " AA' NZ : 2.200e+01\n", - " Factor NZ : 4.500e+01\n", - " Factor Ops : 2.850e+02 (less than 1 second per iteration)\n", + " AA' NZ : 1.700e+01\n", + " Factor NZ : 2.800e+01\n", + " Factor Ops : 1.400e+02 (less than 1 second per iteration)\n", " Threads : 1\n", "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 -8.36906934e-01 -8.36906934e-01 1.57e+00 7.34e-01 6.19e-01 0s\n", - " 1 -9.44994564e-01 -1.48469043e+00 1.36e-01 8.08e-07 8.01e-02 0s\n", - " 2 -7.95781338e-01 -8.91732825e-01 1.50e-07 1.21e-08 8.00e-03 0s\n", - " 3 -8.15171524e-01 -8.29882203e-01 1.65e-13 1.33e-14 1.23e-03 0s\n", - " 4 -8.24009582e-01 -8.24184428e-01 3.18e-14 3.20e-16 1.46e-05 0s\n", - " 5 -8.24036160e-01 -8.24039018e-01 1.54e-13 7.22e-15 2.38e-07 0s\n", - " 6 -8.24036608e-01 -8.24036837e-01 1.10e-10 1.94e-13 1.90e-08 0s\n", + " 0 6.73751896e-01 1.79000000e+00 1.90e+00 6.81e-01 1.00e+00 0s\n", + " 1 1.20057128e+00 2.19855901e+00 1.99e-01 7.50e-07 1.45e-01 0s\n", + " 2 1.11075775e+00 1.33966849e+00 2.19e-07 8.25e-13 2.29e-02 0s\n", + " 3 1.17894103e+00 1.20876338e+00 3.10e-13 2.22e-16 2.98e-03 0s\n", + " 4 1.19359138e+00 1.19386663e+00 1.07e-14 2.22e-16 2.75e-05 0s\n", + " 5 1.19381906e+00 1.19382014e+00 3.32e-12 6.66e-16 1.08e-07 0s\n", "\n", - "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n", - "Optimal objective -8.24036608e-01\n", + "Barrier solved model in 5 iterations and 0.02 seconds (0.00 work units)\n", + "Optimal objective 1.19381906e+00\n", "\n", - "w = [0.50000 0.36395 0.13605],\n", - "z = 0.04810269528692041.\n" + "w = [0.50000 0.32623 0.17377],\n", + "z = 0.82431.\n" ] } ], @@ -500,26 +510,25 @@ "model = gp.Model(\"robustGS\")\n", "\n", "# define variables of interest as a continuous\n", - "w = model.addMVar(shape=problem_size, vtype=GRB.CONTINUOUS, name=\"w\")\n", - "z = model.addVar(name=\"z\")\n", + "w = model.addMVar(shape=problem_size, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "z = model.addVar(lb=0.0, name=\"z\")\n", "\n", "# setup the robust objective function\n", "model.setObjective(\n", - " 0.5*w@(relationship_matrix@w) - lam*w.transpose()@expected_breeding_values - lam*kappa*z,\n", - "GRB.MINIMIZE)\n", + " w.transpose()@expected_breeding_values - (lam/2)*w@(relationship_matrix@w) - kappa*z,\n", + "GRB.MAXIMIZE)\n", "\n", - "# add quadratic uncertainty constraint\n", - "model.addConstr(z**2 <= np.inner(w, omega@w), name=\"uncertainty\")\n", - "model.addConstr(z >= 0, name=\"z positive\")\n", + "# add quadratic uncertainty constraint⁵\n", + "model.addConstr(z**2 >= w.transpose()@omega@w, name=\"uncertainty\")\n", "# add sub-to-half constraints\n", "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", - "# add weight-bound constraints~\n", - "model.addConstr(w >= lower_bound, name=\"lower bound\")\n", - "model.addConstr(w <= upper_bound, name=\"upper bound\")\n", + "# add weight-bound constraints\n", + "model.addConstr(w >= lower_bound, name=\"lower-bound\")\n", + "model.addConstr(w <= upper_bound, name=\"upper-bound\")\n", "\n", "# solve the problem with Gurobi\n", "model.optimize()\n", - "print(f\"w = {w.X},\\nz = {z.X}.\")" + "print(f\"w = {w.X},\\nz = {z.X:.5f}.\")" ] }, { @@ -528,20 +537,23 @@ "source": [ "We can repeat this experiment with varying $\\kappa$ to see how how different tolerances for uncertainty impact the robust turning point returned.\n", "\n", - "| $\\kappa$ | $w$ | $z$ | $f(w,z)$ |\n", - "| ---: | ------------------------: | ------: | -------: |\n", - "| 0.0 | [0.50000 0.37500 0.12500] | 0.02756 | -0.81250 |\n", - "| 0.5 | [0.50000 0.36395 0.13605] | 0.04810 | -0.82404 |\n", - "| 1.0 | [0.50000 0.35282 0.14718] | 0.05204 | -0.83655 |\n", - "| 2.0 | [0.50000 0.33073 0.16927] | 0.05985 | -0.86451 |\n", - "| 4.0 | [0.50000 0.28663 0.21337] | 0.07544 | -0.93214 |\n", - "| 8.0 | [0.50000 0.19816 0.30184] | 0.10671 | -1.11428 |\n", - "| 16.0 | [0.50000 0.07511 0.42489] | 0.15022 | -1.65453 |\n", - "| 32.0 | [0.50000 0.07511 0.42489] | 0.15022 | -2.85630 |\n", - "\n", - "We can see that in the case $\\kappa = 0$ we return the standard optimization solution, as expected. However, we do have $z\\neq0$ which suggests that a small amount of numerical error has been introduced as a result of asking Gurobi to optimize $z$ as well as $w$.\n", - "\n", - "We also note that there are not further changes once $\\kappa$ falls above a certain threshold." + "| $\\kappa$ | $w$ | $z$ | $f(w,z)$ |\n", + "| -------: | ------------------------: | ------: | -------: |\n", + "| 0.00 | [0.50000 0.37500 0.12500] | 3.52875 | 1.62499880e+00\n", + "| 0.25 | [0.50000 0.34998 0.15002] | 0.86184 | 1.40453330e+00\n", + "| 0.50 | [0.50000 0.32623 0.17377] | 0.82431 | 1.19381906e+00\n", + "| 0.75 | [0.50000 0.30443 0.19557] | 0.79089 | 9.91998750e-01\n", + "| 1.00 | [0.50000 0.28391 0.21609] | 0.75044 | 7.98188968e-01\n", + "| 2.00 | [0.50000 0.21877 0.28123] | 0.67181 | 8.61190139e-02\n", + "| 4.00 | [0.50000 0.14644 0.35356] | 0.59279 | -1.16409198e+00\n", + "| 8.00 | [0.50000 0.09102 0.40898] | 0.55140 | -3.43139120e+00\n", + "| 16.00 | [0.50000 0.05658 0.44342] | 0.53608 | -7.76342696e+00\n", + "| 32.00 | [0.50000 0.03695 0.46305] | 0.53128 | -1.62903237e+01\n", + "| 64.00 | [0.50000 0.02635 0.47365] | 0.52992 | -3.32626440e+01\n", + "\n", + "We can see that in the case $\\kappa = 0$ we return the standard optimization solution, as expected. However, we do have $z\\neq0$ which suggests that a small amount of numerical error has been introduced as a result of asking Gurobi to optimize $z$ as well as $w$. We also not that as $\\kappa\\rightarrow\\infty$, we tend towards the portfolio of maximum risk.\n", + "\n", + "TODO _do a comparison of the value $z$ against $\\sqrt{w\\Omega w}$, analyzing the gap between the two. I suspect this will explain away what I originally thought was \"numerical error\"_" ] }, { @@ -648,11 +660,11 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 4 rows, 50 columns and 100 nonzeros\n", - "Model fingerprint: 0x534ff85e\n", + "Model fingerprint: 0x32c74476\n", "Model has 1275 quadratic objective terms\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", - " Objective range [7e-01, 1e+00]\n", + " Objective range [1e+00, 2e+00]\n", " QObjective range [5e-02, 1e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [5e-01, 5e-01]\n", @@ -671,100 +683,97 @@ "\n", " Objective Residual\n", "Iter Primal Dual Primal Dual Compl Time\n", - " 0 3.76397773e+05 -4.22619805e+05 2.50e+04 1.12e+00 9.99e+05 0s\n", - " 1 -5.01117677e+01 -9.99230833e+02 2.68e+01 1.25e-04 1.09e+03 0s\n", - " 2 -8.81632395e-01 -9.63667643e+02 2.68e-05 1.25e-10 1.93e+01 0s\n", - " 3 -8.81606636e-01 -2.03034027e+00 5.19e-09 2.41e-14 2.30e-02 0s\n", - " 4 -9.04332054e-01 -1.04180295e+00 4.73e-10 2.19e-15 2.75e-03 0s\n", - " 5 -9.59656693e-01 -1.01318293e+00 1.67e-15 8.88e-16 1.07e-03 0s\n", - " 6 -9.73793052e-01 -9.76849441e-01 2.22e-16 1.33e-15 6.11e-05 0s\n", - " 7 -9.76163966e-01 -9.76239391e-01 1.39e-15 6.66e-16 1.51e-06 0s\n", - " 8 -9.76213879e-01 -9.76217890e-01 9.91e-15 5.55e-16 8.02e-08 0s\n", - " 9 -9.76215229e-01 -9.76215592e-01 1.65e-13 5.55e-16 7.25e-09 0s\n", - " 10 -9.76215438e-01 -9.76215482e-01 7.03e-13 6.66e-16 8.86e-10 0s\n", - " 11 -9.76215473e-01 -9.76215476e-01 4.24e-13 6.66e-16 4.91e-11 0s\n", + " 0 -3.30175742e+05 4.22619805e+05 2.50e+04 2.24e+00 9.98e+05 0s\n", + " 1 1.02710649e+02 9.99185369e+02 2.74e+01 1.27e-04 1.11e+03 0s\n", + " 2 1.80860414e+00 9.64341109e+02 2.74e-05 1.27e-10 1.93e+01 0s\n", + " 3 1.80859983e+00 3.14837910e+00 1.07e-08 4.98e-14 2.68e-02 0s\n", + " 4 1.87921786e+00 2.10382674e+00 7.18e-10 3.34e-15 4.49e-03 0s\n", + " 5 2.00889156e+00 2.08732673e+00 7.77e-16 1.11e-15 1.57e-03 0s\n", + " 6 2.03435346e+00 2.03715345e+00 3.89e-16 1.78e-15 5.60e-05 0s\n", + " 7 2.03638420e+00 2.03645510e+00 2.22e-15 3.55e-15 1.42e-06 0s\n", + " 8 2.03642992e+00 2.03643695e+00 8.99e-14 2.00e-15 1.41e-07 0s\n", + " 9 2.03643447e+00 2.03643516e+00 2.21e-13 1.55e-15 1.38e-08 0s\n", + " 10 2.03643494e+00 2.03643500e+00 6.20e-13 2.22e-15 1.31e-09 0s\n", + " 11 2.03643499e+00 2.03643499e+00 1.10e-12 1.33e-15 8.34e-11 0s\n", "\n", - "Barrier solved model in 11 iterations and 0.04 seconds (0.00 work units)\n", - "Optimal objective -9.76215473e-01\n", + "Barrier solved model in 11 iterations and 0.02 seconds (0.00 work units)\n", + "Optimal objective 2.03643499e+00\n", "\n", "Gurobi Optimizer version 11.0.0 build v11.0.0rc2 (linux64 - \"Ubuntu 22.04.4 LTS\")\n", "\n", "CPU model: Intel(R) Core(TM) i5-8350U CPU @ 1.70GHz, instruction set [SSE2|AVX|AVX2]\n", "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", - "Optimize a model with 5 rows, 51 columns and 101 nonzeros\n", - "Model fingerprint: 0x8c48581e\n", + "Optimize a model with 4 rows, 51 columns and 100 nonzeros\n", + "Model fingerprint: 0x806a4c54\n", "Model has 1275 quadratic objective terms\n", - "Model has 50 quadratic constraints\n", + "Model has 1 quadratic constraint\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", - " QMatrix range [7e-06, 1e+00]\n", - " Objective range [7e-01, 1e+00]\n", + " QMatrix range [1e-05, 1e+00]\n", + " Objective range [1e+00, 1e+01]\n", " QObjective range [5e-02, 1e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [5e-01, 5e-01]\n", - "Presolve removed 3 rows and 0 columns\n", + "Presolve removed 2 rows and 1 columns\n", "\n", "Continuous model is non-convex -- solving as a MIP\n", "\n", - "Presolve removed 3 rows and 0 columns\n", - "Presolve time: 0.01s\n", - "Presolved: 2552 rows, 1327 columns, 7600 nonzeros\n", + "Presolve removed 2 rows and 1 columns\n", + "Presolve time: 0.00s\n", + "Presolved: 2503 rows, 1325 columns, 6325 nonzeros\n", "Presolved model has 1275 quadratic objective terms\n", - "Presolved model has 1 quadratic constraint(s)\n", "Presolved model has 1275 bilinear constraint(s)\n", - "Variable types: 1327 continuous, 0 integer (0 binary)\n", - "Found heuristic solution: objective -0.9763584\n", + "Variable types: 1325 continuous, 0 integer (0 binary)\n", + "Found heuristic solution: objective 2.0364350\n", "\n", - "Root relaxation: objective -1.285097e+00, 31 iterations, 0.00 seconds (0.00 work units)\n", + "Root relaxation: interrupted, 37 iterations, 0.00 seconds (0.00 work units)\n", "\n", " Nodes | Current Node | Objective Bounds | Work\n", " Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time\n", "\n", - " 0 0 -1.28510 0 36 -0.97636 -1.28510 31.6% - 0s\n", - " 0 0 -1.13066 0 69 -0.97636 -1.13066 15.8% - 0s\n", - " 0 0 -1.05344 0 51 -0.97636 -1.05344 7.89% - 0s\n", - " 0 0 - 0 -0.97636 -1.02501 4.98% - 0s\n", - "\n", - "Cutting planes:\n", - " RLT: 6\n", + " 0 0 - 0 2.03643 2.03643 0.00% - 0s\n", "\n", - "Explored 1 nodes (1828 simplex iterations) in 0.29 seconds (0.10 work units)\n", + "Explored 1 nodes (37 simplex iterations) in 0.03 seconds (0.01 work units)\n", "Thread count was 8 (of 8 available processors)\n", "\n", - "Solution count 1: -0.976358 \n", + "Solution count 1: 2.03643 \n", "\n", "Optimal solution found (tolerance 5.00e-02)\n", - "Best objective -9.763583540500e-01, best bound -1.025012483088e+00, gap 4.9832%\n", + "Best objective 2.036434950287e+00, best bound 2.036434994087e+00, gap 0.0000%\n", "\n", "SIRE WEIGHTS\t\t\t DAM WEIGHTS\n", "--------------------\t\t --------------------\n", " i w_std w_rbs\t\t i w_std w_rbs\n", - "00 0.00000 0.00002 01 0.08319 0.08339\n", - "02 0.11649 0.11508 03 0.00000 0.00001\n", - "04 0.00000 0.00002 05 0.00000 0.00002\n", - "06 0.00000 0.00001 07 0.06817 0.06830\n", - "08 0.00000 0.00001 09 0.03762 0.03701\n", - "10 0.00000 0.00002 11 0.00000 0.00002\n", - "12 0.00000 0.00001 13 0.00355 0.00564\n", - "14 0.11542 0.11450 15 0.06966 0.06913\n", - "16 0.00000 0.00006 17 0.00000 0.00001\n", - "18 0.00000 0.00002 19 0.00000 0.00003\n", - "20 0.00000 0.00002 21 0.00001 0.00062\n", - "22 0.00000 0.00017 23 0.00000 0.00004\n", - "24 0.03537 0.03493 25 0.00000 0.00001\n", - "26 0.00000 0.00001 27 0.00000 0.00002\n", - "28 0.00000 0.00002 29 0.00000 0.00002\n", - "30 0.00000 0.00001 31 0.03580 0.03612\n", - "32 0.08447 0.08517 33 0.00000 0.00003\n", - "34 0.00000 0.00001 35 0.02849 0.02747\n", - "36 0.13489 0.13487 37 0.07446 0.07264\n", - "38 0.00000 0.00001 39 0.00000 0.00001\n", - "40 0.00000 0.00003 41 0.07312 0.07388\n", - "42 0.00000 0.00004 43 0.00000 0.00003\n", - "44 0.01336 0.01485 45 0.00000 0.00002\n", - "46 0.00000 0.00009 47 0.02593 0.02555\n", - "48 0.00000 0.00003 49 0.00000 0.00000\n" + "00 0.00000 0.00000 01 0.10276 0.10276\n", + "02 0.14509 0.14509 03 0.00000 0.00000\n", + "04 0.00000 0.00000 05 0.00000 0.00000\n", + "06 0.00000 0.00000 07 0.10172 0.10172\n", + "08 0.00000 0.00000 09 0.02222 0.02222\n", + "10 0.00000 0.00000 11 0.00000 0.00000\n", + "12 0.00000 0.00000 13 0.00000 0.00000\n", + "14 0.11196 0.11196 15 0.07168 0.07168\n", + "16 0.00000 0.00000 17 0.00000 0.00000\n", + "18 0.00000 0.00000 19 0.00000 0.00000\n", + "20 0.00000 0.00000 21 0.00000 0.00000\n", + "22 0.00000 0.00000 23 0.00000 0.00000\n", + "24 0.00000 0.00000 25 0.00000 0.00000\n", + "26 0.00000 0.00000 27 0.00000 0.00000\n", + "28 0.00000 0.00000 29 0.00000 0.00000\n", + "30 0.00000 0.00000 31 0.02845 0.02845\n", + "32 0.07177 0.07177 33 0.00000 0.00000\n", + "34 0.00000 0.00000 35 0.00000 0.00000\n", + "36 0.17118 0.17118 37 0.08277 0.08277\n", + "38 0.00000 0.00000 39 0.00000 0.00000\n", + "40 0.00000 0.00000 41 0.09038 0.09038\n", + "42 0.00000 0.00000 43 0.00000 0.00000\n", + "44 0.00000 0.00000 45 0.00000 0.00000\n", + "46 0.00000 0.00000 47 0.00001 0.00002\n", + "48 0.00000 0.00000 49 0.00000 0.00000\n", + "\n", + "Maximum change: 0.00001\n", + "Average change: 0.00000\n", + "Minimum change: 0.00000\n" ] } ], @@ -773,10 +782,11 @@ " \"Example/A50.txt\",\n", " \"Example/EBV50.txt\",\n", " \"Example/S50.txt\",\n", - " 50)\n", + " 50\n", + ")\n", "\n", - "lam = 0.5\n", - "kappa = 2\n", + "lam = 1\n", + "kappa = 10\n", "\n", "# define the M so that column i is [1;0] if i is a sire (so even) and [0;1] otherwise \n", "M = np.zeros((2, n))\n", @@ -790,29 +800,28 @@ "model_rbs = gp.Model(\"n50robust\")\n", "\n", "# initialise w for both models, z for robust model\n", - "w_std = model_std.addMVar(shape=n, vtype=GRB.CONTINUOUS, name=\"w\") \n", - "w_rbs = model_rbs.addMVar(shape=n, vtype=GRB.CONTINUOUS, name=\"w\")\n", - "z_rbs = model_rbs.addVar(name=\"z\")\n", + "w_std = model_std.addMVar(shape=n, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\") \n", + "w_rbs = model_rbs.addMVar(shape=n, lb=0.0, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "z_rbs = model_rbs.addVar(lb=0.0, name=\"z\")\n", "\n", "# define the objective functions for both models\n", "model_std.setObjective(\n", - " 0.5*w_std@(sigma@w_std) - lam*w_std.transpose()@mubar,\n", - "GRB.MINIMIZE)\n", + " w_std.transpose() @ mubar - (lam/2)*w_std.transpose()@(sigma@w_std),\n", + "GRB.MAXIMIZE)\n", "\n", "model_rbs.setObjective(\n", " # Gurobi does offer a way to set one objective in terms of another, i.e.\n", " # we could use `model_std.getObjective() - lam*kappa*z_rbs` to define this\n", " # robust objective, but it results in a significant slowdown in code.\n", - " 0.5*w_rbs@(sigma@w_rbs) - lam*w_rbs.transpose()@mubar - lam*kappa*z_rbs,\n", - "GRB.MINIMIZE)\n", + " w_rbs.transpose() @ mubar - (lam/2)*w_rbs.transpose()@(sigma@w_rbs) - kappa*z_rbs,\n", + "GRB.MAXIMIZE)\n", "\n", "# add sum-to-half constraints to both models\n", "model_std.addConstr(M @ w_std == m, name=\"sum-to-half\")\n", - "model_rbs.addConstr(M @ w_std == m, name=\"sum-to-half\")\n", + "model_rbs.addConstr(M @ w_rbs == m, name=\"sum-to-half\")\n", "\n", "# add quadratic uncertainty constraint to the robust model\n", - "model_rbs.addConstr(z_rbs**2 <= np.inner(w_rbs, omega@w_rbs), name=\"uncertainty\")\n", - "model_rbs.addConstr(z_rbs >= 0, name=\"z positive\")\n", + "model_rbs.addConstr(z_rbs**2 <= w_rbs.transpose()@omega@w_rbs, name=\"uncertainty\")\n", "\n", "# since working with non-trivial size, set a time limit\n", "time_limit = 60*5 # 5 minutes\n", @@ -824,6 +833,10 @@ "model_std.setParam('MIPGap', duality_gap)\n", "model_rbs.setParam('MIPGap', duality_gap)\n", "\n", + "# produce model outputs for external usage\n", + "model_std.write(\"MeanVar.mps\")\n", + "model_rbs.write(\"Robust.mps\")\n", + "\n", "# solve both problems with Gurobi\n", "model_std.optimize()\n", "model_rbs.optimize()\n", @@ -834,14 +847,20 @@ "print(\" i w_std w_rbs\\t\\t i w_std w_rbs\")\n", "for candidate in range(25):\n", " print(f\"{candidate*2:02d} {w_std.X[candidate*2]:.5f} {w_rbs.X[candidate*2]:.5f} \\\n", - " {candidate*2+1:02d} {w_std.X[candidate*2+1]:.5f} {w_rbs.X[candidate*2+1]:.5f}\")" + " {candidate*2+1:02d} {w_std.X[candidate*2+1]:.5f} {w_rbs.X[candidate*2+1]:.5f}\")\n", + " \n", + "print(f\"\\nMaximum change: {max(np.abs(w_std.X-w_rbs.X)):.5f}\")\n", + "print(f\"Average change: {np.mean(np.abs(w_std.X-w_rbs.X)):.5f}\")\n", + "print(f\"Minimum change: {min(np.abs(w_std.X-w_rbs.X)):.5f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "We can see that with parameters $\\lambda = 0.5, \\kappa = 2$, we go from having most candidates having zero contribution in the standard solution to most candidates having a contribution of order 0.00001 in the robust solution." + "We can see that with parameters $\\lambda = 1, \\kappa = 10$, we go from having most candidates having zero contribution in the standard solution to most candidates having a contribution of order 0.0001 in the robust solution.\n", + "\n", + "BUG: _with issue regarding Gurobi and `np.inner` fixed, there's now little-if-any difference between the standard and robust solutions. Need to look at what's going on here._" ] }, { @@ -850,19 +869,25 @@ "source": [ "## Footnotes\n", "\n", - "1. In reality, $\\Omega$ is also an estimated variable but we can ignore this for now to avoid going down a rabbit hole of uncertainties.\n", + "1.\n", + " Working in terms of target risk ($\\sigma_p$) rather than target return ($\\mu_p$) means we cannot work through the full range of $\\lambda$ as easily. The \"end point\" of the frontier which previously corresponded to $\\lambda = 0$ now instead corresponds to $\\lambda = \\infty$ (although practically, a value of $\\lambda = 10^5$ recovers it).\n", + "\n", + "2.\n", + " In reality, $\\Omega$ is also an estimated variable but we can ignore this for now to avoid going down a rabbit hole of uncertainties.\n", "\n", - "2. For example, we could alternatively define $U_{\\mu}$ using a box uncertainty set\n", + "3.\n", + " For example, we could alternatively define $U_{\\mu}$ using a box uncertainty set\n", " $$\n", " U_{\\mu} := \\left\\lbrace \\mu :\\ | \\mu_i - \\bar{\\mu}_i| \\leq\\xi_i,\\ \\forall i \\right\\rbrace\n", " $$\n", " where $\\xi$ is some upper bound on the uncertainty value. In this case, work by [Heckel et al (2016)](https://www.risk.net/journal-of-investment-strategies/2475975/insights-into-robust-optimization-decomposing-into-mean-variance-and-risk-based-portfolios) shows that our objective function would gain an additional absolute value term, becoming\n", " $$\n", - " \\min_w \\frac{1}{2}w^{T}\\Sigma w - \\lambda w^{T}\\mu - \\kappa\\|\\Omega^{\\frac{1}{2}} w\\|\\ \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", + " \\max_w \\left(w^{T}\\mu - \\max_i(\\xi_i|w_i|) - \\frac{\\lambda}{2}w^{T}\\Sigma w\\right) \\text{ subject to }\\ Mw = \\begin{bmatrix} 0.5 \\\\ 0.5\\end{bmatrix},\\ l\\leq w\\leq u.\n", " $$\n", " This absolute value term means we lose many of the properties that we would like to work with in our objective functions.\n", "\n", - "3. By \"ball\" we do not just mean a sphere in the typical geometric sense, but a ball in the sense of metric spaces. If we recognise that in our definition we have the transpose of a vector, multiplied by a positive definite matrix, multiplied by the same vector, we can see that\n", + "4.\n", + " By \"ball\" we do not just mean a sphere in the typical geometric sense, but a ball in the sense of metric spaces. If we recognise that in our definition we have the transpose of a vector, multiplied by a positive definite matrix, multiplied by the same vector, we can see that\n", " \\begin{align*}\n", " {(\\mu-\\bar{\\mu})}^{T}\\Omega^{-1}(\\mu-\\bar{\\mu}) &=: x^{T}Ax \\\\\n", " &= x^{T}L^{T}Lx \\quad \\text{(using Cholesky factorisation)} \\\\\n", @@ -873,7 +898,16 @@ " \n", " In other words, in the norm of $\\Omega$'s Cholesky factor our uncertainty set is a ball of radius $\\kappa$ centred on $\\bar{\\mu}$.\n", "\n", - "3. Some fourth footnote which has yet to be written. When is, increment and retain. " + "5.\n", + " Due to issues with Gurobi, it's important that our uncertainty constraint is modelled as\n", + " ```python3\n", + " model.addConstr(z**2 <= w.transpose()@omega@w, name=\"uncertainty\")\n", + " ```\n", + " rather than something like\n", + " ```python3\n", + " model.addConstr(z**2 <= np.inner(w, omega@w), name=\"uncertainty\")\n", + " ```\n", + " For reasons unknown, Gurobi treats these two constraints differently." ] } ],