From 595b7f2f2340c076c4db3cd56f932465e9085e08 Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Tue, 23 Apr 2024 12:06:24 +0100 Subject: [PATCH 1/8] Initial reframing as return maximization --- robust-genetics.ipynb | 316 ++++++++++++++++++++++-------------------- 1 file changed, 164 insertions(+), 152 deletions(-) diff --git a/robust-genetics.ipynb b/robust-genetics.ipynb index 178193f..152118f 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)}\")" ] }, { @@ -240,7 +247,7 @@ "\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" ] @@ -255,30 +262,30 @@ "$$\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,12 +330,12 @@ "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", @@ -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", + "Optimal objective 1.62500000e+00\n", "\n", "w = [0.50000 0.37500 0.12500]\n" ] @@ -371,8 +378,8 @@ "\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", @@ -389,19 +396,19 @@ "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", - " \\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", + " \\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", - "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", + "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", - " \\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", + " \\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", - "TODO: _better explain the idea of $z$ pushing up so that this switch doesn't have any practical difference._\n", + "TODO: _better explain the idea of $z$ pushing down so that this switch doesn't have any practical difference._\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", "$$" ] }, @@ -455,43 +462,44 @@ "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", + "Model fingerprint: 0xd0398316\n", "Model has 3 quadratic objective terms\n", "Model has 3 quadratic constraints\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 time: 0.01s\n", - "Presolved: 9 rows, 8 columns, 14 nonzeros\n", - "Presolved model has 3 second-order cone constraints\n", + "Presolved: 6 rows, 6 columns, 9 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.100e+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 -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.75206688e-01 1.66071429e+00 1.95e+00 7.60e-01 1.40e+00 0s\n", + " 1 9.18871990e-01 2.22754663e+00 2.08e-01 8.36e-07 1.99e-01 0s\n", + " 2 1.26521618e+00 1.51422516e+00 1.85e-02 9.20e-13 3.06e-02 0s\n", + " 3 1.27457950e+00 1.28458784e+00 2.03e-08 4.44e-16 1.11e-03 0s\n", + " 4 1.28096283e+00 1.28138761e+00 8.88e-13 6.66e-16 4.72e-05 0s\n", + " 5 1.28122931e+00 1.28125496e+00 3.34e-12 5.11e-15 2.85e-06 0s\n", + " 6 1.28124988e+00 1.28125699e+00 1.16e-10 2.63e-13 7.90e-07 0s\n", + " 7 1.28125000e+00 1.28125001e+00 1.41e-09 4.38e-12 9.12e-10 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 7 iterations and 0.03 seconds (0.00 work units)\n", + "Optimal objective 1.28125000e+00\n", "\n", - "w = [0.50000 0.36395 0.13605],\n", - "z = 0.04810269528692041.\n" + "w = [0.50000 0.31249 0.18751],\n", + "z = 0.62499.\n" ] } ], @@ -505,11 +513,11 @@ "\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**2 >= np.inner(w, omega@w), name=\"uncertainty\")\n", "model.addConstr(z >= 0, name=\"z positive\")\n", "# add sub-to-half constraints\n", "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", @@ -519,7 +527,7 @@ "\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,19 +536,20 @@ "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", + "| $\\kappa$ | $w$ | $z$ | $f(w,z)$ |\n", + "| -------: | ------------------------: | ------: | -------: |\n", + "| 0.00 | [0.50000 0.37500 0.12500] | 3.47856 | 1.62500 |\n", + "| 0.25 | [0.50000 0.34374 0.15626] | 0.68748 | 1.44531 |\n", + "| 0.50 | [0.50000 0.31249 0.18751] | 0.62499 | 1.28125 |\n", + "| 0.75 | [0.50000 0.28140 0.21860] | 0.56280 | 1.13281 |\n", + "| 1.00 | [0.50000 0.25032 0.24968] | 0.50063 | 1.00000 |\n", + "| 2.00 | [0.50000 0.25000 0.25000] | 0.50000 | 0.50000 |\n", + "| 4.00 | [0.50000 0.25000 0.25000] | 0.50000 | -0.50000 |\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", + "TODO _return to the numerical error notebook and verify whether this is still the case._\n", + "\n", "We also note that there are not further changes once $\\kappa$ falls above a certain threshold." ] }, @@ -648,11 +657,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,21 +680,21 @@ "\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.03 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", @@ -693,13 +702,13 @@ "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", + "Model fingerprint: 0x1b4e6980\n", "Model has 1275 quadratic objective terms\n", "Model has 50 quadratic constraints\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", + " Objective range [1e+00, 8e+00]\n", " QObjective range [5e-02, 1e+00]\n", " Bounds range [0e+00, 0e+00]\n", " RHS range [5e-01, 5e-01]\n", @@ -714,57 +723,55 @@ "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", + "Found heuristic solution: objective 2.0353573\n", "\n", - "Root relaxation: objective -1.285097e+00, 31 iterations, 0.00 seconds (0.00 work units)\n", + "Root relaxation: interrupted, 25 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", + " 0 0 - 0 2.03536 2.03643 0.05% - 0s\n", "\n", - "Cutting planes:\n", - " RLT: 6\n", - "\n", - "Explored 1 nodes (1828 simplex iterations) in 0.29 seconds (0.10 work units)\n", + "Explored 1 nodes (25 simplex iterations) in 0.07 seconds (0.03 work units)\n", "Thread count was 8 (of 8 available processors)\n", "\n", - "Solution count 1: -0.976358 \n", + "Solution count 1: 2.03536 \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.035357260351e+00, best bound 2.036434994087e+00, gap 0.0530%\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.00010 01 0.10276 0.10151\n", + "02 0.14509 0.14458 03 0.00000 0.00000\n", + "04 0.00000 0.00005 05 0.00000 0.00020\n", + "06 0.00000 0.00012 07 0.10172 0.10018\n", + "08 0.00000 0.00017 09 0.02222 0.02115\n", + "10 0.00000 0.00011 11 0.00000 0.00009\n", + "12 0.00000 0.00010 13 0.00000 0.00059\n", + "14 0.11196 0.11075 15 0.07168 0.07044\n", + "16 0.00000 0.00004 17 0.00000 0.00006\n", + "18 0.00000 0.00008 19 0.00000 0.00010\n", + "20 0.00000 0.00010 21 0.00000 0.00038\n", + "22 0.00000 0.00011 23 0.00000 0.00019\n", + "24 0.00000 0.00093 25 0.00000 0.00023\n", + "26 0.00000 0.00011 27 0.00000 0.00012\n", + "28 0.00000 0.00006 29 0.00000 0.00014\n", + "30 0.00000 0.00011 31 0.02845 0.02736\n", + "32 0.07177 0.07176 33 0.00000 0.00039\n", + "34 0.00000 0.00012 35 0.00000 0.00335\n", + "36 0.17118 0.16967 37 0.08277 0.07995\n", + "38 0.00000 0.00007 39 0.00000 0.00016\n", + "40 0.00000 0.00010 41 0.09038 0.08835\n", + "42 0.00000 0.00009 43 0.00000 0.00018\n", + "44 0.00000 0.00035 45 0.00000 0.00008\n", + "46 0.00000 0.00021 47 0.00001 0.00474\n", + "48 0.00000 0.00011 49 0.00000 0.00005\n", + "\n", + "Maximum change: 0.00473\n", + "Average change: 0.00057\n", + "Minimum change: 0.00000\n" ] } ], @@ -773,10 +780,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", @@ -796,19 +804,19 @@ "\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,\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", @@ -834,14 +842,18 @@ "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." ] }, { @@ -873,7 +885,7 @@ " \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. " + "4. Note that working in terms of target risk $\\sigma_p$ rather than target return $\\mu_p$ means we can no longer work through the full range of $\\lambda$, since the end of the \"end\" of the critical frontier now corresponds to the point when $\\lambda = \\infty$ (though practically, a value of $\\lambda = 10^5$ recovers it). " ] } ], From 217a157220df42a0ab2d6295848d421ae4395c17 Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Tue, 23 Apr 2024 13:31:37 +0100 Subject: [PATCH 2/8] Updated footnotes Includes a fix where in one place I'd written in the norm-form of a term coming from the quadratic uncertainty set, rather than the term coming from the box uncertainty set --- robust-genetics.ipynb | 32 ++++++++++++++++++++------------ 1 file changed, 20 insertions(+), 12 deletions(-) diff --git a/robust-genetics.ipynb b/robust-genetics.ipynb index 152118f..6d12f1a 100644 --- a/robust-genetics.ipynb +++ b/robust-genetics.ipynb @@ -212,7 +212,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "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.⁴" + "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.¹" ] }, { @@ -243,7 +243,7 @@ "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", @@ -258,11 +258,11 @@ "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 $\\min_{\\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", " \\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", @@ -398,7 +398,11 @@ "source": [ "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", + " \\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", + "or in its norm form\n", + "$$\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", "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", @@ -862,19 +866,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. For example, we could alternatively define $U_{\\mu}$ using a box uncertainty set\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", + "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", @@ -883,9 +893,7 @@ " &= {\\|x\\|}_{L^{-1}}^2 \\leq \\kappa^2. \n", " \\end{align*}\n", " \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", - "4. Note that working in terms of target risk $\\sigma_p$ rather than target return $\\mu_p$ means we can no longer work through the full range of $\\lambda$, since the end of the \"end\" of the critical frontier now corresponds to the point when $\\lambda = \\infty$ (though practically, a value of $\\lambda = 10^5$ recovers it). " + " In other words, in the norm of $\\Omega$'s Cholesky factor our uncertainty set is a ball of radius $\\kappa$ centred on $\\bar{\\mu}$." ] } ], From a44afc213c82427854f891702e2873222589a0c5 Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Tue, 23 Apr 2024 13:56:15 +0100 Subject: [PATCH 3/8] Added notebook with alternate framing Has the original working for risk minimization mindset, rather than the return maximization mindset --- Scratch/risk-minimization.ipynb | 819 ++++++++++++++++++++++++++++++++ 1 file changed, 819 insertions(+) create mode 100644 Scratch/risk-minimization.ipynb diff --git a/Scratch/risk-minimization.ipynb b/Scratch/risk-minimization.ipynb new file mode 100644 index 0000000..293d2ae --- /dev/null +++ b/Scratch/risk-minimization.ipynb @@ -0,0 +1,819 @@ +{ + "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, 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", + "TODO: _better explain the idea of $z$ pushing up so that this switch doesn't have any practical difference._\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 23 rows, 4 columns and 25 nonzeros\n", + "Model fingerprint: 0xb99edb8a\n", + "Model has 3 quadratic objective terms\n", + "Model has 3 quadratic constraints\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 22 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", + "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", + " 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", + "\n", + "Barrier solved model in 6 iterations and 0.02 seconds (0.00 work units)\n", + "Optimal objective -8.24036608e-01\n", + "\n", + "w = [0.50000 0.36395 0.13605],\n", + "z = 0.04810269528692041.\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, vtype=GRB.CONTINUOUS, name=\"w\")\n", + "z = model.addVar(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 <= np.inner(w, omega@w), name=\"uncertainty\")\n", + "model.addConstr(z >= 0, name=\"z positive\")\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)$ |\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." + ] + }, + { + "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": 11, + "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.03 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 5 rows, 51 columns and 101 nonzeros\n", + "Model fingerprint: 0x8c48581e\n", + "Model has 1275 quadratic objective terms\n", + "Model has 50 quadratic constraints\n", + "Coefficient statistics:\n", + " Matrix range [1e+00, 1e+00]\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + " QMatrix range [7e-06, 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 3 rows and 0 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", + "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", + "\n", + "Root relaxation: objective -1.285097e+00, 31 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", + "\n", + "Explored 1 nodes (1828 simplex iterations) in 0.26 seconds (0.10 work units)\n", + "Thread count was 8 (of 8 available processors)\n", + "\n", + "Solution count 1: -0.976358 \n", + "\n", + "Optimal solution found (tolerance 5.00e-02)\n", + "Best objective -9.763583540500e-01, best bound -1.025012483088e+00, gap 4.9832%\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", + "\n", + "Maximum change: 0.00209\n", + "Average change: 0.00029\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, 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", + "\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 <= np.inner(w_rbs, omega@w_rbs), name=\"uncertainty\")\n", + "model_rbs.addConstr(z_rbs >= 0, name=\"z positive\")\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$, 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." + ] + } + ], + "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 +} From d4bacccacc93d51c5dd34e84f8ae3fa31e1eda83 Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Tue, 23 Apr 2024 14:19:06 +0100 Subject: [PATCH 4/8] Updated Julian's $n = 2$E --- Example2/n2-test.ipynb | 123 ++++++++++++++++++++--------------------- 1 file changed, 61 insertions(+), 62 deletions(-) diff --git a/Example2/n2-test.ipynb b/Example2/n2-test.ipynb index c4fb745..d649094 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", + "Barrier solved model in 7 iterations and 0.03 seconds (0.00 work units)\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", @@ -197,43 +197,42 @@ "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", + "Model fingerprint: 0x41322026\n", "Model has 2 quadratic objective terms\n", "Model has 2 quadratic constraints\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", + "Presolve time: 0.04s\n", + "Presolved: 7 rows, 8 columns, 12 nonzeros\n", "Presolved model has 3 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.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 -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 -2.28881626e-02 1.22656250e+00 9.39e-01 1.43e-01 3.26e-01 0s\n", + " 1 6.91883948e-01 1.06699842e+00 1.98e-01 1.57e-07 5.84e-02 0s\n", + " 2 6.44324078e-01 7.21385315e-01 2.17e-07 1.73e-13 7.71e-03 0s\n", + " 3 6.67762956e-01 6.69433547e-01 6.21e-09 6.66e-16 1.67e-04 0s\n", + " 4 6.68355745e-01 6.68371340e-01 1.53e-13 2.00e-15 1.56e-06 0s\n", + " 5 6.68366724e-01 6.68367382e-01 3.57e-11 6.05e-13 6.58e-08 0s\n", "\n", - "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n", - "Optimal objective -6.73611064e-01\n", + "Barrier solved model in 5 iterations and 0.05 seconds (0.00 work units)\n", + "Optimal objective 6.68366724e-01\n", "\n", - "Standard: [0.25000 0.75000]\n", - "Robust: [0.41660 0.58340]\n" + "Standard: [0.00000001 0.99999999]\n", + "Robust: [0.85714285 0.14285715]\n" ] } ], @@ -249,13 +248,13 @@ "z_rbs = model_rbs.addVar(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**2 >= np.inner(w_rbs, omega@w_rbs), name=\"uncertainty\")\n", "model_rbs.addConstr(z_rbs >= 0, name=\"z positive\")\n", "\n", "# solve problem with Gurobi, print alongside standard solution for comparison\n", @@ -268,7 +267,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 eight decimal places." ] }, { @@ -286,43 +285,43 @@ "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", + "Model fingerprint: 0xe31c34b4\n", "Model has 2 quadratic objective terms\n", "Model has 2 quadratic constraints\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", + "Presolve time: 0.02s\n", + "Presolved: 7 rows, 8 columns, 12 nonzeros\n", "Presolved model has 3 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.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 -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.04e+00 6.81e-01 5.08e-01 0s\n", + " 1 1.54200502e+00 1.52822126e+00 4.89e-02 1.42e-01 7.24e-02 0s\n", + " 2 1.64909037e+00 1.63852570e+00 5.38e-08 3.21e-02 1.66e-02 0s\n", + " 3 1.74359514e+00 1.75057927e+00 3.23e-10 8.27e-05 7.48e-04 0s\n", + " 4 1.74999141e+00 1.75002072e+00 1.01e-12 9.47e-08 2.99e-06 0s\n", + " 5 1.74999468e+00 1.75000778e+00 3.11e-12 3.00e-08 1.33e-06 0s\n", + " 6 1.74999882e+00 1.75000092e+00 1.06e-11 1.95e-09 2.11e-07 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 6 iterations and 0.06 seconds (0.00 work units)\n", + "Optimal objective 1.74999882e+00\n", "\n", - "Standard: [0.2500000031 0.7499999969]\n", - "Robust: [0.2499847190 0.7500152810]\n" + "Standard: [0.0000000076 0.9999999924]\n", + "Robust: [0.0000000230 0.9999999770]\n" ] } ], @@ -338,13 +337,13 @@ "z_rbs = model_rbs.addVar(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**2 >= np.inner(w_rbs, omega@w_rbs), name=\"uncertainty\")\n", "model_rbs.addConstr(z_rbs >= 0, name=\"z positive\")\n", "\n", "# solve problem with Gurobi, print alongside standard solution for comparison\n", From f0bc21dee251bccf9d29266c7ba3fead4282b381 Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Tue, 23 Apr 2024 15:58:24 +0100 Subject: [PATCH 5/8] Added a better explanation of $z$ pushing downwards And similarly pushing upwards in the risk minimization framing --- Scratch/risk-minimization.ipynb | 2 +- robust-genetics.ipynb | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/Scratch/risk-minimization.ipynb b/Scratch/risk-minimization.ipynb index 293d2ae..d26083f 100644 --- a/Scratch/risk-minimization.ipynb +++ b/Scratch/risk-minimization.ipynb @@ -344,7 +344,7 @@ "$$\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", - "TODO: _better explain the idea of $z$ pushing up so that this switch doesn't have any practical difference._\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", diff --git a/robust-genetics.ipynb b/robust-genetics.ipynb index 6d12f1a..97cb310 100644 --- a/robust-genetics.ipynb +++ b/robust-genetics.ipynb @@ -408,7 +408,7 @@ "$$\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", - "TODO: _better explain the idea of $z$ pushing down so that this switch doesn't have any practical difference._\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\\geq\\sqrt{w^{T}\\Omega w}$ can be squared on both sides:\n", "$$\n", From d70735f5af8cd89e29bc21ed5affa6a1f43d0a4a Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Tue, 23 Apr 2024 18:08:58 +0100 Subject: [PATCH 6/8] Removed a TODO I'd todone --- robust-genetics.ipynb | 6 ++---- 1 file changed, 2 insertions(+), 4 deletions(-) diff --git a/robust-genetics.ipynb b/robust-genetics.ipynb index 97cb310..7eb5310 100644 --- a/robust-genetics.ipynb +++ b/robust-genetics.ipynb @@ -550,11 +550,9 @@ "| 2.00 | [0.50000 0.25000 0.25000] | 0.50000 | 0.50000 |\n", "| 4.00 | [0.50000 0.25000 0.25000] | 0.50000 | -0.50000 |\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", + "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 note that there are not further changes once $\\kappa$ falls above a certain threshold.\n", "\n", - "TODO _return to the numerical error notebook and verify whether this is still the case._\n", - "\n", - "We also note that there are not further changes once $\\kappa$ falls above a certain threshold." + "TODO _do a comparison of the value $z$ against $\\sqrt{w\\Omega w}$, analyzing the gap between the two._" ] }, { From c854716301078c1d468b0b57631cb44c27f3af9a Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Thu, 25 Apr 2024 11:44:58 +0100 Subject: [PATCH 7/8] Incorporate bounds into variable definitions --- Example2/n2-test.ipynb | 28 +++++------ Scratch/numerical-error.ipynb | 25 +++++----- Scratch/risk-minimization.ipynb | 42 +++++++---------- robust-genetics.ipynb | 82 +++++++++++++++++---------------- 4 files changed, 83 insertions(+), 94 deletions(-) diff --git a/Example2/n2-test.ipynb b/Example2/n2-test.ipynb index d649094..b93e15f 100644 --- a/Example2/n2-test.ipynb +++ b/Example2/n2-test.ipynb @@ -196,8 +196,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 2 rows, 3 columns and 3 nonzeros\n", - "Model fingerprint: 0x41322026\n", + "Optimize a model with 1 rows, 3 columns and 2 nonzeros\n", + "Model fingerprint: 0xd120428b\n", "Model has 2 quadratic objective terms\n", "Model has 2 quadratic constraints\n", "Coefficient statistics:\n", @@ -207,8 +207,7 @@ " 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.04s\n", + "Presolve time: 0.01s\n", "Presolved: 7 rows, 8 columns, 12 nonzeros\n", "Presolved model has 3 second-order cone constraints\n", "Ordering time: 0.00s\n", @@ -228,7 +227,7 @@ " 4 6.68355745e-01 6.68371340e-01 1.53e-13 2.00e-15 1.56e-06 0s\n", " 5 6.68366724e-01 6.68367382e-01 3.57e-11 6.05e-13 6.58e-08 0s\n", "\n", - "Barrier solved model in 5 iterations and 0.05 seconds (0.00 work units)\n", + "Barrier solved model in 5 iterations and 0.03 seconds (0.00 work units)\n", "Optimal objective 6.68366724e-01\n", "\n", "Standard: [0.00000001 0.99999999]\n", @@ -244,8 +243,8 @@ "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", " w_rbs.transpose()@mubar - (lam/2)*w_rbs.transpose()@(sigma@w_rbs) - kappa*z_rbs,\n", @@ -255,7 +254,6 @@ "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", "\n", "# solve problem with Gurobi, print alongside standard solution for comparison\n", "model_rbs.optimize()\n", @@ -284,8 +282,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 2 rows, 3 columns and 3 nonzeros\n", - "Model fingerprint: 0xe31c34b4\n", + "Optimize a model with 1 rows, 3 columns and 2 nonzeros\n", + "Model fingerprint: 0x6c4430f5\n", "Model has 2 quadratic objective terms\n", "Model has 2 quadratic constraints\n", "Coefficient statistics:\n", @@ -295,8 +293,7 @@ " 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.02s\n", + "Presolve time: 0.01s\n", "Presolved: 7 rows, 8 columns, 12 nonzeros\n", "Presolved model has 3 second-order cone constraints\n", "Ordering time: 0.00s\n", @@ -317,7 +314,7 @@ " 5 1.74999468e+00 1.75000778e+00 3.11e-12 3.00e-08 1.33e-06 0s\n", " 6 1.74999882e+00 1.75000092e+00 1.06e-11 1.95e-09 2.11e-07 0s\n", "\n", - "Barrier solved model in 6 iterations and 0.06 seconds (0.00 work units)\n", + "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n", "Optimal objective 1.74999882e+00\n", "\n", "Standard: [0.0000000076 0.9999999924]\n", @@ -333,8 +330,8 @@ "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", " w_rbs.transpose()@mubar - (lam/2)*w_rbs.transpose()@(sigma@w_rbs) - kappa*z_rbs,\n", @@ -344,7 +341,6 @@ "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", "\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..791a356 100644 --- a/Scratch/numerical-error.ipynb +++ b/Scratch/numerical-error.ipynb @@ -150,7 +150,7 @@ " 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.03 seconds (0.00 work units)\n", + "Barrier solved model in 8 iterations and 0.04 seconds (0.00 work units)\n", "Optimal objective -8.12500000e-01\n", "\n", "w = [0.50000 0.37500 0.12500]\n" @@ -162,7 +162,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 +172,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 +231,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 +242,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 +263,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 +277,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 +287,11 @@ "\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 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", diff --git a/Scratch/risk-minimization.ipynb b/Scratch/risk-minimization.ipynb index d26083f..337a592 100644 --- a/Scratch/risk-minimization.ipynb +++ b/Scratch/risk-minimization.ipynb @@ -302,7 +302,7 @@ " 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", + "Barrier solved model in 8 iterations and 0.03 seconds (0.00 work units)\n", "Optimal objective -8.12500000e-01\n", "\n", "w = [0.50000 0.37500 0.12500]\n" @@ -314,7 +314,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", @@ -401,8 +401,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: 0xb99edb8a\n", + "Optimize a model with 22 rows, 4 columns and 24 nonzeros\n", + "Model fingerprint: 0x5b10d4aa\n", "Model has 3 quadratic objective terms\n", "Model has 3 quadratic constraints\n", "Coefficient statistics:\n", @@ -412,7 +412,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", @@ -447,8 +447,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", @@ -457,7 +457,6 @@ "\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 sub-to-half constraints\n", "model.addConstr(M @ w == m, name=\"sum-to-half\")\n", "# add weight-bound constraints\n", @@ -568,7 +567,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 10, "metadata": {}, "outputs": [ { @@ -628,28 +627,22 @@ "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: 0x74debc94\n", "Model has 1275 quadratic objective terms\n", "Model has 50 quadratic constraints\n", "Coefficient statistics:\n", - " Matrix range [1e+00, 1e+00]\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ + " Matrix range [1e+00, 1e+00]\n", " QMatrix range [7e-06, 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 3 rows and 0 columns\n", + "Presolve removed 2 rows and 0 columns\n", "\n", "Continuous model is non-convex -- solving as a MIP\n", "\n", - "Presolve removed 3 rows and 0 columns\n", + "Presolve removed 2 rows and 0 columns\n", "Presolve time: 0.01s\n", "Presolved: 2552 rows, 1327 columns, 7600 nonzeros\n", "Presolved model has 1275 quadratic objective terms\n", @@ -671,7 +664,7 @@ "Cutting planes:\n", " RLT: 6\n", "\n", - "Explored 1 nodes (1828 simplex iterations) in 0.26 seconds (0.10 work units)\n", + "Explored 1 nodes (1828 simplex iterations) in 0.22 seconds (0.10 work units)\n", "Thread count was 8 (of 8 available processors)\n", "\n", "Solution count 1: -0.976358 \n", @@ -736,9 +729,9 @@ "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", @@ -758,7 +751,6 @@ "\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", "\n", "# since working with non-trivial size, set a time limit\n", "time_limit = 60*5 # 5 minutes\n", diff --git a/robust-genetics.ipynb b/robust-genetics.ipynb index 7eb5310..d22db74 100644 --- a/robust-genetics.ipynb +++ b/robust-genetics.ipynb @@ -362,7 +362,7 @@ " 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", + "Barrier solved model in 8 iterations and 0.02 seconds (0.00 work units)\n", "Optimal objective 1.62500000e+00\n", "\n", "w = [0.50000 0.37500 0.12500]\n" @@ -374,7 +374,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", @@ -384,8 +384,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", @@ -465,8 +465,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: 0xd0398316\n", + "Optimize a model with 22 rows, 4 columns and 24 nonzeros\n", + "Model fingerprint: 0x5341cd5b\n", "Model has 3 quadratic objective terms\n", "Model has 3 quadratic constraints\n", "Coefficient statistics:\n", @@ -476,7 +476,7 @@ " 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: 6 rows, 6 columns, 9 nonzeros\n", "Presolved model has 2 second-order cone constraints\n", @@ -499,7 +499,7 @@ " 6 1.28124988e+00 1.28125699e+00 1.16e-10 2.63e-13 7.90e-07 0s\n", " 7 1.28125000e+00 1.28125001e+00 1.41e-09 4.38e-12 9.12e-10 0s\n", "\n", - "Barrier solved model in 7 iterations and 0.03 seconds (0.00 work units)\n", + "Barrier solved model in 7 iterations and 0.02 seconds (0.00 work units)\n", "Optimal objective 1.28125000e+00\n", "\n", "w = [0.50000 0.31249 0.18751],\n", @@ -512,8 +512,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", @@ -522,12 +522,11 @@ "\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 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", @@ -703,77 +702,77 @@ "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: 0x1b4e6980\n", + "Optimize a model with 4 rows, 51 columns and 100 nonzeros\n", + "Model fingerprint: 0xacf3991a\n", "Model has 1275 quadratic objective terms\n", "Model has 50 quadratic constraints\n", "Coefficient statistics:\n", " Matrix range [1e+00, 1e+00]\n", " QMatrix range [7e-06, 1e+00]\n", - " Objective range [1e+00, 8e+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 0 columns\n", "\n", "Continuous model is non-convex -- solving as a MIP\n", "\n", - "Presolve removed 3 rows and 0 columns\n", + "Presolve removed 2 rows and 0 columns\n", "Presolve time: 0.01s\n", "Presolved: 2552 rows, 1327 columns, 7600 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 2.0353573\n", + "Found heuristic solution: objective 2.0353486\n", "\n", "Root relaxation: interrupted, 25 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 - 0 2.03536 2.03643 0.05% - 0s\n", + " 0 0 - 0 2.03535 2.03643 0.05% - 0s\n", "\n", "Explored 1 nodes (25 simplex iterations) in 0.07 seconds (0.03 work units)\n", "Thread count was 8 (of 8 available processors)\n", "\n", - "Solution count 1: 2.03536 \n", + "Solution count 1: 2.03535 \n", "\n", "Optimal solution found (tolerance 5.00e-02)\n", - "Best objective 2.035357260351e+00, best bound 2.036434994087e+00, gap 0.0530%\n", + "Best objective 2.035348603347e+00, best bound 2.036434994087e+00, gap 0.0534%\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.00010 01 0.10276 0.10151\n", - "02 0.14509 0.14458 03 0.00000 0.00000\n", + "00 0.00000 0.00010 01 0.10276 0.10148\n", + "02 0.14509 0.14457 03 0.00000 0.00012\n", "04 0.00000 0.00005 05 0.00000 0.00020\n", - "06 0.00000 0.00012 07 0.10172 0.10018\n", + "06 0.00000 0.00012 07 0.10172 0.10016\n", "08 0.00000 0.00017 09 0.02222 0.02115\n", "10 0.00000 0.00011 11 0.00000 0.00009\n", "12 0.00000 0.00010 13 0.00000 0.00059\n", - "14 0.11196 0.11075 15 0.07168 0.07044\n", + "14 0.11196 0.11075 15 0.07168 0.07042\n", "16 0.00000 0.00004 17 0.00000 0.00006\n", "18 0.00000 0.00008 19 0.00000 0.00010\n", - "20 0.00000 0.00010 21 0.00000 0.00038\n", + "20 0.00000 0.00009 21 0.00000 0.00038\n", "22 0.00000 0.00011 23 0.00000 0.00019\n", "24 0.00000 0.00093 25 0.00000 0.00023\n", "26 0.00000 0.00011 27 0.00000 0.00012\n", "28 0.00000 0.00006 29 0.00000 0.00014\n", - "30 0.00000 0.00011 31 0.02845 0.02736\n", - "32 0.07177 0.07176 33 0.00000 0.00039\n", + "30 0.00000 0.00011 31 0.02845 0.02737\n", + "32 0.07177 0.07175 33 0.00000 0.00039\n", "34 0.00000 0.00012 35 0.00000 0.00335\n", - "36 0.17118 0.16967 37 0.08277 0.07995\n", + "36 0.17118 0.16970 37 0.08277 0.07995\n", "38 0.00000 0.00007 39 0.00000 0.00016\n", - "40 0.00000 0.00010 41 0.09038 0.08835\n", + "40 0.00000 0.00010 41 0.09038 0.08832\n", "42 0.00000 0.00009 43 0.00000 0.00018\n", "44 0.00000 0.00035 45 0.00000 0.00008\n", - "46 0.00000 0.00021 47 0.00001 0.00474\n", + "46 0.00000 0.00020 47 0.00001 0.00472\n", "48 0.00000 0.00011 49 0.00000 0.00005\n", "\n", - "Maximum change: 0.00473\n", - "Average change: 0.00057\n", - "Minimum change: 0.00000\n" + "Maximum change: 0.00471\n", + "Average change: 0.00058\n", + "Minimum change: 0.00002\n" ] } ], @@ -800,9 +799,9 @@ "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", @@ -813,7 +812,7 @@ " # 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", - " w_rbs.transpose() @ mubar - (lam/2)*w_rbs.transpose()@(sigma@w_rbs) - kappa*z,\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", @@ -822,7 +821,6 @@ "\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", "\n", "# since working with non-trivial size, set a time limit\n", "time_limit = 60*5 # 5 minutes\n", @@ -834,6 +832,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", From c82f19e5df5e80cd03408a80b7a11d47a243b2d5 Mon Sep 17 00:00:00 2001 From: Joshua Fogg Date: Thu, 25 Apr 2024 13:15:49 +0100 Subject: [PATCH 8/8] Bypassed Gurobi's issue with `np.inner` --- Example2/n2-test.ipynb | 74 ++++++------- Scratch/numerical-error.ipynb | 16 ++- Scratch/risk-minimization.ipynb | 191 +++++++++++++++++--------------- robust-genetics.ipynb | 176 +++++++++++++++-------------- 4 files changed, 247 insertions(+), 210 deletions(-) diff --git a/Example2/n2-test.ipynb b/Example2/n2-test.ipynb index b93e15f..345b6cd 100644 --- a/Example2/n2-test.ipynb +++ b/Example2/n2-test.ipynb @@ -142,7 +142,7 @@ " 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.03 seconds (0.00 work units)\n", + "Barrier solved model in 7 iterations and 0.02 seconds (0.00 work units)\n", "Optimal objective 1.75000000e+00\n", "\n", "Standard: [0.00000001 0.99999999]\n" @@ -197,9 +197,9 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 1 rows, 3 columns and 2 nonzeros\n", - "Model fingerprint: 0xd120428b\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", @@ -208,30 +208,31 @@ " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+00, 1e+00]\n", "Presolve time: 0.01s\n", - "Presolved: 7 rows, 8 columns, 12 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.700e+01\n", - " Factor NZ : 2.800e+01\n", - " Factor Ops : 1.400e+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 -2.28881626e-02 1.22656250e+00 9.39e-01 1.43e-01 3.26e-01 0s\n", - " 1 6.91883948e-01 1.06699842e+00 1.98e-01 1.57e-07 5.84e-02 0s\n", - " 2 6.44324078e-01 7.21385315e-01 2.17e-07 1.73e-13 7.71e-03 0s\n", - " 3 6.67762956e-01 6.69433547e-01 6.21e-09 6.66e-16 1.67e-04 0s\n", - " 4 6.68355745e-01 6.68371340e-01 1.53e-13 2.00e-15 1.56e-06 0s\n", - " 5 6.68366724e-01 6.68367382e-01 3.57e-11 6.05e-13 6.58e-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 5 iterations and 0.03 seconds (0.00 work units)\n", - "Optimal objective 6.68366724e-01\n", + "Barrier solved model in 6 iterations and 0.03 seconds (0.00 work units)\n", + "Optimal objective 5.52216328e-01\n", "\n", "Standard: [0.00000001 0.99999999]\n", - "Robust: [0.85714285 0.14285715]\n" + "Robust: [0.83155134 0.16844866]\n" ] } ], @@ -253,7 +254,7 @@ "# 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**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", @@ -265,7 +266,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note that if $\\kappa = 0$, the standard solution not _exactly_ recovered, only accurate to eight decimal places." + "Note that if $\\kappa = 0$, the standard solution not _exactly_ recovered, only accurate to seven decimal places." ] }, { @@ -283,9 +284,9 @@ "Thread count: 4 physical cores, 8 logical processors, using up to 8 threads\n", "\n", "Optimize a model with 1 rows, 3 columns and 2 nonzeros\n", - "Model fingerprint: 0x6c4430f5\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", @@ -294,31 +295,30 @@ " Bounds range [0e+00, 0e+00]\n", " RHS range [1e+00, 1e+00]\n", "Presolve time: 0.01s\n", - "Presolved: 7 rows, 8 columns, 12 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.700e+01\n", - " Factor NZ : 2.800e+01\n", - " Factor Ops : 1.400e+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 1.22656250e+00 1.22656250e+00 1.04e+00 6.81e-01 5.08e-01 0s\n", - " 1 1.54200502e+00 1.52822126e+00 4.89e-02 1.42e-01 7.24e-02 0s\n", - " 2 1.64909037e+00 1.63852570e+00 5.38e-08 3.21e-02 1.66e-02 0s\n", - " 3 1.74359514e+00 1.75057927e+00 3.23e-10 8.27e-05 7.48e-04 0s\n", - " 4 1.74999141e+00 1.75002072e+00 1.01e-12 9.47e-08 2.99e-06 0s\n", - " 5 1.74999468e+00 1.75000778e+00 3.11e-12 3.00e-08 1.33e-06 0s\n", - " 6 1.74999882e+00 1.75000092e+00 1.06e-11 1.95e-09 2.11e-07 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 1.74999882e+00\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.0000000076 0.9999999924]\n", - "Robust: [0.0000000230 0.9999999770]\n" + "Robust: [0.0000001531 0.9999998469]\n" ] } ], @@ -340,7 +340,7 @@ "# 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**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 791a356..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." ] }, { @@ -150,7 +159,7 @@ " 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.04 seconds (0.00 work units)\n", + "Barrier solved model in 8 iterations and 0.03 seconds (0.00 work units)\n", "Optimal objective -8.12500000e-01\n", "\n", "w = [0.50000 0.37500 0.12500]\n" @@ -287,6 +296,7 @@ "\n", "# add quadratic uncertainty constraint\n", "model.addConstr(z**2 <= np.inner(w, omega@w), name=\"uncertainty\")\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", @@ -302,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 index 337a592..7a6d8c0 100644 --- a/Scratch/risk-minimization.ipynb +++ b/Scratch/risk-minimization.ipynb @@ -302,7 +302,7 @@ " 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.03 seconds (0.00 work units)\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" @@ -402,9 +402,9 @@ "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: 0x5b10d4aa\n", + "Model fingerprint: 0xf7819a35\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", @@ -413,32 +413,46 @@ " 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: 9 rows, 8 columns, 14 nonzeros\n", - "Presolved model has 3 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", - " Threads : 1\n", + "Continuous model is non-convex -- solving as a MIP\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", - "\n", - "Barrier solved model in 6 iterations and 0.02 seconds (0.00 work units)\n", - "Optimal objective -8.24036608e-01\n", - "\n", - "w = [0.50000 0.36395 0.13605],\n", - "z = 0.04810269528692041.\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" ] } ], @@ -456,7 +470,7 @@ "GRB.MINIMIZE)\n", "\n", "# add quadratic uncertainty constraint\n", - "model.addConstr(z**2 <= np.inner(w, omega@w), name=\"uncertainty\")\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", @@ -474,20 +488,15 @@ "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", + "| $\\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 also note that there are not further changes once $\\kappa$ falls above a certain threshold." + "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." ] }, { @@ -619,7 +628,7 @@ " 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.03 seconds (0.00 work units)\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", @@ -628,12 +637,12 @@ "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: 0x74debc94\n", + "Model fingerprint: 0xd47c065e\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", + " 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", @@ -643,66 +652,70 @@ "Continuous model is non-convex -- solving as a MIP\n", "\n", "Presolve removed 2 rows and 0 columns\n", - "Presolve time: 0.01s\n", - "Presolved: 2552 rows, 1327 columns, 7600 nonzeros\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 -0.9763584\n", + "Found heuristic solution: objective -1.1617295\n", "\n", - "Root relaxation: objective -1.285097e+00, 31 iterations, 0.00 seconds (0.00 work units)\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 -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", + " 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", - "Cutting planes:\n", - " RLT: 6\n", - "\n", - "Explored 1 nodes (1828 simplex iterations) in 0.22 seconds (0.10 work units)\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: -0.976358 \n", + "Solution count 1: -1.16173 \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 -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.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", - "\n", - "Maximum change: 0.00209\n", - "Average change: 0.00029\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" ] } @@ -750,7 +763,7 @@ "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**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", @@ -783,7 +796,7 @@ "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 = 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." ] } ], diff --git a/robust-genetics.ipynb b/robust-genetics.ipynb index d22db74..1507c41 100644 --- a/robust-genetics.ipynb +++ b/robust-genetics.ipynb @@ -339,7 +339,7 @@ " 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", @@ -362,7 +362,7 @@ " 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.02 seconds (0.00 work units)\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" @@ -466,9 +466,9 @@ "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: 0x5341cd5b\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", @@ -478,32 +478,30 @@ " RHS range [5e-01, 1e+00]\n", "Presolve removed 21 rows and 1 columns\n", "Presolve time: 0.01s\n", - "Presolved: 6 rows, 6 columns, 9 nonzeros\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 : 1.100e+01\n", - " Factor NZ : 2.100e+01\n", - " Factor Ops : 9.100e+01 (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 6.75206688e-01 1.66071429e+00 1.95e+00 7.60e-01 1.40e+00 0s\n", - " 1 9.18871990e-01 2.22754663e+00 2.08e-01 8.36e-07 1.99e-01 0s\n", - " 2 1.26521618e+00 1.51422516e+00 1.85e-02 9.20e-13 3.06e-02 0s\n", - " 3 1.27457950e+00 1.28458784e+00 2.03e-08 4.44e-16 1.11e-03 0s\n", - " 4 1.28096283e+00 1.28138761e+00 8.88e-13 6.66e-16 4.72e-05 0s\n", - " 5 1.28122931e+00 1.28125496e+00 3.34e-12 5.11e-15 2.85e-06 0s\n", - " 6 1.28124988e+00 1.28125699e+00 1.16e-10 2.63e-13 7.90e-07 0s\n", - " 7 1.28125000e+00 1.28125001e+00 1.41e-09 4.38e-12 9.12e-10 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 7 iterations and 0.02 seconds (0.00 work units)\n", - "Optimal objective 1.28125000e+00\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.31249 0.18751],\n", - "z = 0.62499.\n" + "w = [0.50000 0.32623 0.17377],\n", + "z = 0.82431.\n" ] } ], @@ -520,11 +518,11 @@ " 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", + "# 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", + "# add weight-bound constraints\n", "model.addConstr(w >= lower_bound, name=\"lower-bound\")\n", "model.addConstr(w <= upper_bound, name=\"upper-bound\")\n", "\n", @@ -541,17 +539,21 @@ "\n", "| $\\kappa$ | $w$ | $z$ | $f(w,z)$ |\n", "| -------: | ------------------------: | ------: | -------: |\n", - "| 0.00 | [0.50000 0.37500 0.12500] | 3.47856 | 1.62500 |\n", - "| 0.25 | [0.50000 0.34374 0.15626] | 0.68748 | 1.44531 |\n", - "| 0.50 | [0.50000 0.31249 0.18751] | 0.62499 | 1.28125 |\n", - "| 0.75 | [0.50000 0.28140 0.21860] | 0.56280 | 1.13281 |\n", - "| 1.00 | [0.50000 0.25032 0.24968] | 0.50063 | 1.00000 |\n", - "| 2.00 | [0.50000 0.25000 0.25000] | 0.50000 | 0.50000 |\n", - "| 4.00 | [0.50000 0.25000 0.25000] | 0.50000 | -0.50000 |\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 note that there are not further changes once $\\kappa$ falls above a certain threshold.\n", - "\n", - "TODO _do a comparison of the value $z$ against $\\sqrt{w\\Omega w}$, analyzing the gap between the two._" + "| 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\"_" ] }, { @@ -694,7 +696,7 @@ " 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.03 seconds (0.00 work units)\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", @@ -703,76 +705,75 @@ "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: 0xacf3991a\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", + " 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 2 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 2 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 2.0353486\n", + "Variable types: 1325 continuous, 0 integer (0 binary)\n", + "Found heuristic solution: objective 2.0364350\n", "\n", - "Root relaxation: interrupted, 25 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 - 0 2.03535 2.03643 0.05% - 0s\n", + " 0 0 - 0 2.03643 2.03643 0.00% - 0s\n", "\n", - "Explored 1 nodes (25 simplex iterations) in 0.07 seconds (0.03 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: 2.03535 \n", + "Solution count 1: 2.03643 \n", "\n", "Optimal solution found (tolerance 5.00e-02)\n", - "Best objective 2.035348603347e+00, best bound 2.036434994087e+00, gap 0.0534%\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.00010 01 0.10276 0.10148\n", - "02 0.14509 0.14457 03 0.00000 0.00012\n", - "04 0.00000 0.00005 05 0.00000 0.00020\n", - "06 0.00000 0.00012 07 0.10172 0.10016\n", - "08 0.00000 0.00017 09 0.02222 0.02115\n", - "10 0.00000 0.00011 11 0.00000 0.00009\n", - "12 0.00000 0.00010 13 0.00000 0.00059\n", - "14 0.11196 0.11075 15 0.07168 0.07042\n", - "16 0.00000 0.00004 17 0.00000 0.00006\n", - "18 0.00000 0.00008 19 0.00000 0.00010\n", - "20 0.00000 0.00009 21 0.00000 0.00038\n", - "22 0.00000 0.00011 23 0.00000 0.00019\n", - "24 0.00000 0.00093 25 0.00000 0.00023\n", - "26 0.00000 0.00011 27 0.00000 0.00012\n", - "28 0.00000 0.00006 29 0.00000 0.00014\n", - "30 0.00000 0.00011 31 0.02845 0.02737\n", - "32 0.07177 0.07175 33 0.00000 0.00039\n", - "34 0.00000 0.00012 35 0.00000 0.00335\n", - "36 0.17118 0.16970 37 0.08277 0.07995\n", - "38 0.00000 0.00007 39 0.00000 0.00016\n", - "40 0.00000 0.00010 41 0.09038 0.08832\n", - "42 0.00000 0.00009 43 0.00000 0.00018\n", - "44 0.00000 0.00035 45 0.00000 0.00008\n", - "46 0.00000 0.00020 47 0.00001 0.00472\n", - "48 0.00000 0.00011 49 0.00000 0.00005\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.00471\n", - "Average change: 0.00058\n", - "Minimum change: 0.00002\n" + "Maximum change: 0.00001\n", + "Average change: 0.00000\n", + "Minimum change: 0.00000\n" ] } ], @@ -820,7 +821,7 @@ "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**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", @@ -857,7 +858,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "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." + "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._" ] }, { @@ -893,7 +896,18 @@ " &= {\\|x\\|}_{L^{-1}}^2 \\leq \\kappa^2. \n", " \\end{align*}\n", " \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}$." + " 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", + "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." ] } ],