{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Risk measures\n",
"\n",
"The value-at-risk (VaR) of a random variable $X$ at risk level $\\alpha$ is defined as\n",
"$$\\textrm{VaR}_\\alpha(X) = \\min_\\eta\\left(\\mathbb{P}[X\\geq\\eta]\\leq 1-\\alpha\\right)$$\n",
"and the conditional value-at-risk (CVaR) is the conditional expected value\n",
"$$\\textrm{CVaR}_\\alpha(X)=\\mathbb{E}(X~|~X\\geq\\textrm{VaR}_\\alpha(X)).$$\n",
"In the standard setting where the random variable $X$ represents a loss (on an investment) and we take $\\alpha=0.9$, CVaR can be interpreted as the expected loss assuming we end up in the $10\\%$ worst scenarios. The concept was studied in the paper Optimization of Conditional Value-at-Risk, Rockafellar and Uryasev, Journal of Risk, 2002, where we also find an equivalent formula\n",
"$$\\textrm{CVaR}_\\alpha(X)=\\min_\\eta\\left(\\eta+(1-\\alpha)^{-1}\\mathbb{E}([X-\\eta]^+)\\right)$$\n",
"where $[x]^+=\\max(x,0)$.\n",
"\n",
"Some authors consider generalizations\n",
"$$\\rho^f_\\alpha(X)=\\min_\\eta\\left(\\eta+(1-\\alpha)^{-1}f^{-1}\\left(\\mathbb{E}(f(X-\\eta))\\right)\\right),$$\n",
"where $f$ is any increasing, convex function on $\\mathbb{R}_+$ which is identically zero on $\\mathbb{R}_-$. Examples include:\n",
"* $f(x)=[x]^+$, recovering the previous definition of CVaR,\n",
"* $f(x)=\\exp([t]^+)-1$, log-exponential convex risk measure (LogExpCR), (Vinel, Krokhmal 2017),\n",
"* $f(x)=([x]^+)^q$, higher moment coherent risk measures (HMCR) (Krokhmal 2007).\n",
"\n",
"# Stochastic portfolio optimization\n",
"\n",
"In this notebook we implement the models described in Section 5 of (Vinel, Krokhmal 2017) with various risk measures. We consider the problem of investing in $n$ stocks given $m$ historical return vectors $r_1,\\ldots,r_m\\in\\mathbb{R}^n$ taken with probabilities $p_1,\\ldots,p_m$. We consider the following conditions:\n",
"* total budget constraint:\n",
" $$1^Tx\\leq 1$$\n",
"* minimal expected return:\n",
" $$x^T\\mathbb{E}(\\mu)\\geq \\bar{r}$$\n",
" where $\\mu$ is the uncertain vector of returns, $\\mathbb{E}(\\mu)=\\sum_j p_jr_j$ is its expected value averaged over the historical scenarios, and $\\bar{r}$ is the minimal required expected return level.\n",
"* no short-selling:\n",
" $$x\\geq 0$$\n",
"* the objective represents minimization of CVaR $\\rho^f_\\alpha(-\\mu^Tx)$ of the loss $-\\mu^Tx$, that is\n",
" $$\\textrm{minimize}_{\\eta,x}\\left(\\eta+(1-\\alpha)^{-1}f^{-1}(\\sum_{j=1}^mp_jf(-r_j^Tx-\\eta))\\right)$$\n",
" \n",
"We used a preprocessed file with closing prices from NYSE for 1675 assets over 2050 trading days. Each historical scenario is defined by computing returns of all assets over some period of 10 trading days They are returned in a martix $r$ with $n$ rows and $m$ columns."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"(1675, 2040)\n"
]
}
],
"source": [
"%matplotlib inline\n",
"import matplotlib\n",
"import matplotlib.pyplot as plt\n",
"from mosek.fusion import *\n",
"import numpy as np \n",
"import sys, os, csv\n",
"\n",
"allPrices = np.array(list(csv.reader(open(\"allPrices.txt\"), lineterminator='\\n', quoting=csv.QUOTE_NONNUMERIC)))\n",
"n,m = allPrices.shape\n",
"p1 = allPrices[:,0:m-10]\n",
"p2 = allPrices[:,10:m]\n",
"r = (p2-p1)/p1\n",
"print(r.shape)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# CVaR\n",
"\n",
"We first implement the linear CVaR model with $f(x)=[x]^+$. We can then write the problem as\n",
"\n",
"$$\n",
"\\begin{array}{ll}\n",
"\\mbox{minimize} & \\eta + (1-\\alpha)^{-1} p^Tw\\\\\n",
"\\mbox{subject to} & 1^Tx\\leq 1, \\\\\n",
" & x^Trp\\geq \\bar{r}, \\\\\n",
" & w_j\\geq \\max\\{0, -r_j^Tx-\\eta \\}, \\quad j=1,\\ldots,m,\n",
"\\end{array}\n",
"$$\n",
"\n",
"with variables $x\\in\\mathbb{R}^n, w\\in\\mathbb{R}^m, \\eta\\in\\mathbb{R}$. Below is the Fusion implementation."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"def modelCVaR(alpha, rbar, p, r, m, n): \n",
" M = Model(\"CVaR\")\n",
"\n",
" # Variables\n",
" X = M.variable('X', n, Domain.greaterThan(0.0))\n",
" W = M.variable('W', m, Domain.greaterThan(0.0))\n",
" eta = M.variable('eta')\n",
"\n",
" # The bounds w_j + r_j^Tx + eta \\geq 0\n",
" M.constraint(Expr.add([W, Expr.mul(np.transpose(r), X), Var.vrepeat(eta, m)]), Domain.greaterThan(0.0))\n",
" # Minimal risk\n",
" M.constraint(Expr.dot(X, np.matmul(r,p)), Domain.greaterThan(rbar)) \n",
" # Budget\n",
" M.constraint(Expr.sum(X), Domain.lessThan(1.0))\n",
"\n",
" # Set the objective \n",
" M.objective(ObjectiveSense.Minimize, Expr.add(eta, Expr.mul(1/(1-alpha), Expr.dot(p, W))))\n",
"\n",
" return M"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# LogExpCR\n",
"\n",
"Next we have $f(x)=\\exp([x]^+)-1$. Since $\\sum_j p_j=1$, the model can equivalently be formulated as\n",
"\n",
"$$\n",
"\\begin{array}{ll}\n",
"\\mbox{minimize} & \\eta + (1-\\alpha)^{-1} t\\\\\n",
"\\mbox{subject to} & 1^Tx\\leq 1, \\\\\n",
" & x^Trp\\geq \\bar{r}, \\\\\n",
" & w_j\\geq \\max\\{0, -r_j^Tx-\\eta \\}, \\quad j=1,\\ldots,m, \\\\\n",
" & t\\geq \\log\\left(\\sum_j p_j\\exp(w_j)\\right).\n",
"\\end{array}\n",
"$$\n",
"\n",
"The last bound (log-sum-exp constraint) can be expressed using an auxiliary variable $\\xi\\in\\mathbb{R}^m$:\n",
"\n",
"$$\n",
"\\begin{array}{l}\n",
"p^T\\xi \\leq 1, \\\\\n",
"\\xi_j\\geq \\exp(w_j-t),\n",
"\\end{array}\n",
"$$\n",
"\n",
"where the last inequality is conic-representable as $(\\xi_j,1,w_j-t)\\in K_\\mathrm{exp}$."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"def modelLogExpCR(alpha, rbar, p, r, m, n): \n",
" M = Model(\"LogExpCR\")\n",
"\n",
" # Variables\n",
" X = M.variable('X', n, Domain.greaterThan(0.0))\n",
" W = M.variable('W', m, Domain.greaterThan(0.0))\n",
" eta = M.variable('eta')\n",
" xi = M.variable('xi', m) \n",
" t = M.variable('t') \n",
"\n",
" # The log-sum-exp constraint\n",
" M.constraint(Expr.dot(p, xi), Domain.lessThan(1.0)) \n",
" # Stack allcones together - every row has the form (xi_j, 1, w_j-t)\n",
" # The matrix notation means that every row belongs to the cone\n",
" M.constraint(Expr.hstack(xi, Expr.constTerm(m, 1.0), Expr.sub(W, Var.vrepeat(t, m))), Domain.inPExpCone())\n",
"\n",
" # Continue with the linear constraints\n",
" M.constraint(Expr.add([W, Expr.mul(np.transpose(r), X), Var.vrepeat(eta, m)]), Domain.greaterThan(0.0))\n",
" M.constraint(Expr.dot(X, np.matmul(r,p)), Domain.greaterThan(rbar)) \n",
" M.constraint(Expr.sum(X), Domain.lessThan(1.0))\n",
"\n",
" # Set the objective \n",
" M.objective(ObjectiveSense.Minimize, Expr.add(eta, Expr.mul(1/(1-alpha), t)))\n",
"\n",
" return M"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# HMCR\n",
"\n",
"The model using $f(x)=([x]^+)^q$ is the same as LogExpCR, except that the last constraint becomes\n",
"\n",
"$$t\\geq\\left(\\sum_j p_jw_j^q\\right)^{1/q}.$$\n",
"\n",
"This $q$-norm cone constraint is equivalent to\n",
"\n",
"$$\n",
"\\begin{array}{l}\n",
"t = 1^T\\xi,\\\\\n",
"t^{1-1/q}\\xi_j^{1/q}\\geq |p_j^{1/q}w_j|\n",
"\\end{array}\n",
"$$\n",
"\n",
"where the last row is precisely a conic constraint involving the $3$-dimenional power cone with parameters $(1-1/q,1/q)$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"def modelHMCR(alpha, rbar, p, r, m, n, q): \n",
" M = Model(\"HMCR\")\n",
"\n",
" # Variables\n",
" X = M.variable('X', n, Domain.greaterThan(0.0))\n",
" W = M.variable('W', m, Domain.greaterThan(0.0))\n",
" eta = M.variable('eta')\n",
" xi = M.variable('xi', m) \n",
" t = M.variable('t') \n",
"\n",
" # The power cones representing the q-norm constraint\n",
" M.constraint(Expr.sub(t, Expr.sum(xi)), Domain.equalsTo(0.0)) \n",
" M.constraint(Expr.hstack(Var.vrepeat(t, m), xi, Expr.mulElm(np.power(p, 1.0/q), W)), Domain.inPPowerCone(1.0-1.0/q))\n",
"\n",
" # Continue with the linear constraints\n",
" M.constraint(Expr.add([W, Expr.mul(np.transpose(r), X), Var.vrepeat(eta, m)]), Domain.greaterThan(0.0))\n",
" M.constraint(Expr.dot(X, np.matmul(r,p)), Domain.greaterThan(rbar)) \n",
" M.constraint(Expr.sum(X), Domain.lessThan(1.0))\n",
"\n",
" # Set the objective \n",
" M.objective(ObjectiveSense.Minimize, Expr.add(eta, Expr.mul(1/(1-alpha), t)))\n",
"\n",
" return M"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Example\n",
"\n",
"We solve a few problems with $n=100$ stocks and $m=2000$ sample historical scenarios."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" CVaR: status=ProblemStatus.PrimalAndDualFeasible time=0.20s\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFBZJREFUeJzt3X+QXWddx/H3121aFhjYlkbHbBsSxhDoCBJcCogi8isFnTbWMqSjYxlxOox0RNE4zeCIFmcKxlFxpoN0oIqMtEDpxAyimUqL/iMlG4KUtATSCu1uigTToCM7NAlf/7hn4eZ6N3vu3l+753m/Znb23nOfe+7z3Gfv55z7nGfPicxEklSGHxp3BSRJo2PoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpy3rgr0Oniiy/OTZs2jbsakrSmHDx48FuZuX65cqsu9Ddt2sTs7Oy4qyFJa0pEfL1OOYd3JKkghr4kFcTQl6SCGPqSVBBDX5IKsupm70gl2Htonj37j3Ds5AIbpibZtX0rO7ZNj7taKoChL43Y3kPz7L7rfhZOnQFg/uQCu++6H8Dg19A5vCON2J79R74f+IsWTp1hz/4jY6qRSmLoSyN27ORCT8ulQXJ4R0PjuHV3G6Ymme8S8BumJsdQG5XGPX0NxeK49fzJBZIfjFvvPTQ/7qqN3a7tW5lcN3HWssl1E+zavnVMNVJJDH0NhePWS9uxbZqbr34e01OTBDA9NcnNVz/Pb0EaiVrDOxFxBfBeYAL4QGa+u+PxtwO/DpwGjgO/lplfrx67Dvj9qugfZ+aHBlR3rWKOW5/bjm3ThrzGYtk9/YiYAG4BXgdcBlwbEZd1FDsEzGTm84E7gT+pnnsR8E7gxcDlwDsj4sLBVV+r1VLj045bS+NVZ3jncuBoZj6cmU8AdwBXtRfIzHsz8zvV3c8Cl1S3twN3Z+aJzHwcuBu4YjBV12rmuLW0OtUZ3pkGHm27P0drz30pbwb+8RzP9TttARaHLpy9I60udUI/uizLrgUjfgWYAX62l+dGxPXA9QAbN26sUSWtBY5bS6tPneGdOeDStvuXAMc6C0XEq4F3AFdm5nd7eW5m3pqZM5k5s379slf7kiStUJ3QPwBsiYjNEXE+sBPY114gIrYB76cV+N9se2g/8NqIuLA6gPvaapkkaQyWHd7JzNMRcQOtsJ4AbsvMwxFxEzCbmfuAPcBTgY9HBMAjmXllZp6IiHfR2nAA3JSZJ4bSEknSsiKz6/D82MzMzKQXRpek3kTEwcycWa6c/5ErSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQWqFfkRcERFHIuJoRNzY5fGXR8TnI+J0RFzT8diZiPhC9bNvUBWXJPXuvOUKRMQEcAvwGmAOOBAR+zLzgbZijwBvAn63yyoWMvMFA6irJKlPy4Y+cDlwNDMfBoiIO4CrgO+HfmZ+rXrse0OooyRpQOoM70wDj7bdn6uW1fWkiJiNiM9GxI6eaidJGqg6e/rRZVn28BobM/NYRDwLuCci7s/Mh856gYjrgesBNm7c2MOqJUm9qLOnPwdc2nb/EuBY3RfIzGPV74eBzwDbupS5NTNnMnNm/fr1dVctSepRndA/AGyJiM0RcT6wE6g1CyciLoyIC6rbFwMvo+1YgCRptJYN/cw8DdwA7AceBD6WmYcj4qaIuBIgIl4UEXPAG4D3R8Th6unPBWYj4t+Be4F3d8z6kSSNUGT2Mjw/fDMzMzk7OzvuakjSmhIRBzNzZrly/keuJBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSpIrdCPiCsi4khEHI2IG7s8/vKI+HxEnI6Iazoeuy4ivlr9XDeoikuSerds6EfEBHAL8DrgMuDaiLiso9gjwJuAj3Q89yLgncCLgcuBd0bEhf1XW5K0EnX29C8Hjmbmw5n5BHAHcFV7gcz8WmZ+Efhex3O3A3dn5onMfBy4G7hiAPWWJK1AndCfBh5tuz9XLaujn+dKkgasTuhHl2VZc/21nhsR10fEbETMHj9+vOaqJUm9Oq9GmTng0rb7lwDHaq5/DnhFx3M/01koM28FbgWYmZmpu0GRpFr2Hppnz/4jHDu5wIapSXZt38qObWUOOtTZ0z8AbImIzRFxPrAT2Fdz/fuB10bEhdUB3NdWyyRpJPYemmf3Xfczf3KBBOZPLrD7rvvZe2h+3FUbi2VDPzNPAzfQCusHgY9l5uGIuCkirgSIiBdFxBzwBuD9EXG4eu4J4F20NhwHgJuqZZI0Env2H2Hh1Jmzli2cOsOe/UfGVKPxqjO8Q2Z+CvhUx7I/aLt9gNbQTbfn3gbc1kcdJWnFjp1c6Gl50/kfuZIabcPUZE/Lm87Ql9Rou7ZvZXLdxFnLJtdNsGv71jHVaLxqDe9I0lq1OEvH2Tsthr6kxtuxbbrYkO/k8I4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klSQWqEfEVdExJGIOBoRN3Z5/IKI+Gj1+H0RsalavikiFiLiC9XPXw22+pKkXpy3XIGImABuAV4DzAEHImJfZj7QVuzNwOOZ+WMRsRN4D/DG6rGHMvMFA663JGkF6uzpXw4czcyHM/MJ4A7gqo4yVwEfqm7fCbwqImJw1ZQkDUKd0J8GHm27P1ct61omM08D3waeUT22OSIORcS/RMTPdHuBiLg+ImYjYvb48eM9NUCSVF+d0O+2x541yzwGbMzMbcDbgY9ExNP+X8HMWzNzJjNn1q9fX6NKkqSVqBP6c8ClbfcvAY4tVSYizgOeDpzIzO9m5n8BZOZB4CHg2f1WWpK0MnVC/wCwJSI2R8T5wE5gX0eZfcB11e1rgHsyMyNifXUgmIh4FrAFeHgwVZck9WrZ2TuZeToibgD2AxPAbZl5OCJuAmYzcx/wQeDDEXEUOEFrwwDwcuCmiDgNnAHekpknhtEQSdLyIrNzeH68ZmZmcnZ2dtzVaIy9h+bZs/8Ix04usGFqkl3bt7JjW+dxeElrXUQczMyZ5cotu6evtWvvoXl233U/C6fOADB/coHdd90PsKaD3w2ZtHKehqHB9uw/8v3AX7Rw6gx79h8ZU436t7ghmz+5QPKDDdneQ/Pjrpq0Jhj6DXbs5EJPy9eCJm7IpFEy9Btsw9RkT8vXgiZuyKRRMvQbbNf2rUyumzhr2eS6CXZt3zqmGvWviRsyaZQM/QbbsW2am69+HtNTkwQwPTXJzVc/b00f9GzihkwaJWfvNNyObdNrOuQ7LbbF2TvSyhj6WnOatiGTRsnhHUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCuI8fa0JgzqdsqdlVukMfa16g7ouQFOvLyD1wuEdrXqDOp2yp2WW3NPXAAx7yGRQp1P2tMySe/rq0yiuZDWo0yl7WmbJ0FefRjFkMqjTKXtaZsnhHfVpFEMmgzqdchNPy+xsJPXK0FdfNkxNMt8l4Ac9ZDKo0ymP87TMgw5oZyNpJRzeUV8cMqlnGMc+nI2klTD01ZcmXpJxGIYR0M5G0ko4vKO+eSWr5Q0joEc1tKZmcU9fGoFhTBd1aE0rYeifw95D87zs3few+cZ/4GXvvmegc89VlmEEtENrWgmHd5bgzAgN0rCmizZpNpJGw9BfwrkOvPmHrZVo0rGPQe4UufEYLYd3luDMCGlpg5qNNIrTeOhshv4SPE+LtLRB7RT5vwaj5/DOEnZt33rW11dwZkQdflUfn1G+94OaLuo36tFzT38JzozonV/Vx2fU7/2gZiMN4hu1s+x6457+OTTpwNsoePB7fEb93g9qNlK/36hX6yy71fyN19DXwNT9qj6sD8Rq/qAN2ziGSQaxU9TvxmM17mis1g3RIkNfA1NnnHdYH4jV/kEbtrV8SoZ+Nh6r8ZjAatwQtXNMXwNTZ5x3WLM1Sp8FUuopGVbjLLvVuCFqVyv0I+KKiDgSEUcj4sYuj18QER+tHr8vIja1Pba7Wn4kIrYPruq984DPcNU5+D2sD8Rq/6ANW6kTD1bjxm41bojaRWaeu0DEBPAV4DXAHHAAuDYzH2gr8xvA8zPzLRGxE/jFzHxjRFwG3A5cDmwA/hl4dmae6XydRTMzMzk7O9tzQ9rHc58+uY4IOPmdU9+//fh3ThFAe2sX7091KX/yO6fOGl9cbv3dnvtzz1nPvV8+Xvs5w7i9VBsGWb8647CLr91tCGJRnX7od7399Gm/fb1UO/qpR511DutvcVDvZWcbBtUvw6rrcu/lSrKm189TNxFxMDNnli1XI/RfCvxhZm6v7u8GyMyb28rsr8r8W0ScB3wDWA/c2F62vdxSr7eS0O8czx2kyXUT/NJPTvOJg/NDWf8ojKINk+smltyzHET/dFv/Stc7zj7tbMcw3pthfh66vfYg3svFNgBD69NR9/ti0HduAOo41+dpyderGfp1hnemgUfb7s9Vy7qWyczTwLeBZ9R8bt+6jecOysKpM9x+36NrNvBhNG041/j5IPqn2/pXut5x9mlnO4bx3gzz89DttQfxXi62YZh9Oup+T2AioufAh+Eej6ozeye6LOtsx1Jl6jyXiLgeuB5g48aNNap0tmGP255Z5tvQWjCKNgx7XL1zPf2sd5x92l7vYbw3oz6OMaj3st9616nHqPu9n9cbVj/W2dOfAy5tu38JcGypMtXwztOBEzWfS2bempkzmTmzfv36+rWvDPsAyUR023atLaNoQ68HsHqtU+d6+lnvOPu0vd6D+tsdxjrrGtR7uWFqsq+6r8Z+7+f1htWPdUL/ALAlIjZHxPnATmBfR5l9wHXV7WuAe7J1sGAfsLOa3bMZ2AJ8bjBV/4FuR/CXEh2/lzO5boJrX3xp7fWvRqNow7lmTCw1w6KXOnVb/0rXO84+7WxHL3+7o1xnL689iPdysQ0rrftq7Pd+Xm+YM5CWHd7JzNMRcQOwH5gAbsvMwxFxEzCbmfuADwIfjoijtPbwd1bPPRwRHwMeAE4Dbz3XzJ2V6vyvvkEdwW8vP/PMi9b87J32Noxy9s65/uuyzvu61PpXst5++nQYs3fq/u32s85RzN7p573s1oZeP2uj6veV9PVKsmOY/02+7OydUVvplE1JKtkgZ+9IkhrC0Jekghj6klQQQ1+SCmLoS1JBVt3snYg4Dny9j1VcDHxrQNVZK0psM5TZ7hLbDGW2u9c2PzMzl/3v1lUX+v2KiNk605aapMQ2Q5ntLrHNUGa7h9Vmh3ckqSCGviQVpImhf+u4KzAGJbYZymx3iW2GMts9lDY3bkxfkrS0Ju7pS5KW0JjQX+7i7U0REZdGxL0R8WBEHI6It1XLL4qIuyPiq9XvC8dd10GLiImIOBQRn6zub46I+6o2f7Q69XejRMRURNwZEV+u+vylTe/riPjt6m/7SxFxe0Q8qYl9HRG3RcQ3I+JLbcu69m20/GWVb1+MiBeu9HUbEfrVxdtvAV4HXAZcW12UvYlOA7+Tmc8FXgK8tWrrjcCnM3ML8OnqftO8DXiw7f57gD+v2vw48Oax1Gq43gv8U2Y+B/gJWu1vbF9HxDTwm8BMZv44rdO576SZff03wBUdy5bq29fRuh7JFlpXGXzfSl+0EaEPXA4czcyHM/MJ4A7gqjHXaSgy87HM/Hx1+39ohcA0rfZ+qCr2IWDHeGo4HBFxCfDzwAeq+wG8ErizKtLENj8NeDmt61WQmU9k5kka3te0rvMxWV2F78nAYzSwrzPzX2ldf6TdUn17FfC32fJZYCoifnQlr9uU0B/JBdhXm4jYBGwD7gN+JDMfg9aGAfjh8dVsKP4C+D3ge9X9ZwAnM/N0db+Jff4s4Djw19Ww1gci4ik0uK8zcx74U+ARWmH/beAgze/rRUv17cAyrimhX+sC7E0SEU8FPgH8Vmb+97jrM0wR8QvANzPzYPviLkWb1ufnAS8E3peZ24D/pUFDOd1UY9hXAZuBDcBTaA1tdGpaXy9nYH/vTQn9Whdgb4qIWEcr8P8uM++qFv/n4te96vc3x1W/IXgZcGVEfI3W0N0rae35T1VDANDMPp8D5jLzvur+nbQ2Ak3u61cD/5GZxzPzFHAX8FM0v68XLdW3A8u4poR+nYu3N0I1lv1B4MHM/LO2h9ovTn8d8PejrtuwZObuzLwkMzfR6tt7MvOXgXuBa6pijWozQGZ+A3g0IhavkP0qWtebbmxf0xrWeUlEPLn6W19sc6P7us1SfbsP+NVqFs9LgG8vDgP1LDMb8QO8HvgK8BDwjnHXZ4jt/GlaX+u+CHyh+nk9rTHuTwNfrX5fNO66Dqn9rwA+Wd1+FvA54CjwceCCcddvCO19ATBb9fde4MKm9zXwR8CXgS8BHwYuaGJfA7fTOm5xitae/JuX6ltawzu3VPl2P63ZTSt6Xf8jV5IK0pThHUlSDYa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kF+T9asd7lM1a2egAAAABJRU5ErkJggg==\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"LogExpCR: status=ProblemStatus.PrimalAndDualFeasible time=0.69s\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAFCpJREFUeJzt3X9sXWd9x/H3FzctBgRuqTctbkPCFgLVYARMgbGxjR9NYFObsaKmG1qRKlVoVDDYMjVCGlv5o7BMAyZVjAq6ARoUKFUWMTaro932zyhxCKOkJRA6oHZgBFJ302rRJHz3xz2GW/c6Pte+P+zzvF+SlXvOfc45z3Of649PnvPccyMzkSSV4QnDroAkaXAMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBzhl2BRa78MILc/PmzcOuhiStK4cOHfpBZo4vV27Nhf7mzZuZnp4edjUkaV2JiG/XKefwjiQVxNCXpIIY+pJUEENfkgpi6EtSQdbc7B2pBPsPz7Jv6ijH5+bZODbKnh3b2LV9YtjVUgEMfWnA9h+eZe8d9zJ/6gwAs3Pz7L3jXgCDX33n8I40YPumjv4k8BfMnzrDvqmjQ6qRSmLoSwN2fG6+q/VSLxn60oBtHBvtar3US4a+NGB7dmxjdMPIY9aNbhhhz45tQ6qRSuKFXGnAFi7WOntHw2Doq2+clri0XdsnfC00FIa++sJpidLa5Ji++sJpidLaZOirL5yWKK1NtUI/InZGxNGIOBYRN3R4/u0RcV9EfCUiPh8Rz2h77pqI+Eb1c00vK6+1y2mJ0tq0bOhHxAhwM/Aa4BLg6oi4ZFGxw8BkZj4PuB34i2rbC4B3Ai8GLgXeGRHn9676WqucliitTXXO9C8FjmXmA5n5KHAbcEV7gcy8OzMfqRa/AFxUPd4B3JmZJzPzIeBOYGdvqq61bNf2CW563XOZGBslgImxUW563XO9iCsNWZ3ZOxPAg23LM7TO3JdyLfBPZ9nW3/pCOC1RWnvqhH50WJcdC0a8AZgEfq2bbSPiOuA6gE2bNtWokiRpJeoM78wAF7ctXwQcX1woIl4FvAO4PDN/1M22mXlLZk5m5uT4+HjdukuSulQn9A8CWyNiS0ScC+wGDrQXiIjtwAdpBf73256aAi6LiPOrC7iXVeskSUOw7PBOZp6OiOtphfUIcGtmHomIG4HpzDwA7AOeAnw6IgC+k5mXZ+bJiHgXrT8cADdm5sm+tESStKzI7Dg8PzSTk5M5PT097GpI0roSEYcyc3K5cn4iV5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVpFboR8TOiDgaEcci4oYOz788Ir4UEacj4spFz52JiC9XPwd6VXFJUvfOWa5ARIwANwOvBmaAgxFxIDPvayv2HeCNwB932MV8Zj6/B3WVJK3SsqEPXAocy8wHACLiNuAK4Cehn5nfqp77cR/qKEnqkTrDOxPAg23LM9W6up4YEdMR8YWI2NWpQERcV5WZPnHiRBe7liR1o07oR4d12cUxNmXmJPC7wPsi4ucft7PMWzJzMjMnx8fHu9i1JKkbdUJ/Bri4bfki4HjdA2Tm8erfB4B/BbZ3UT9JUg/VCf2DwNaI2BIR5wK7gVqzcCLi/Ig4r3p8IfAy2q4FSJIGa9nQz8zTwPXAFHA/8KnMPBIRN0bE5QAR8aKImAFeD3wwIo5Umz8HmI6I/wTuBt69aNaPJGmAIrOb4fn+m5yczOnp6WFXQ5LWlYg4VF0/PSs/kStJBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klSQWqEfETsj4mhEHIuIGzo8//KI+FJEnI6IKxc9d01EfKP6uaZXFZckdW/Z0I+IEeBm4DXAJcDVEXHJomLfAd4IfHzRthcA7wReDFwKvDMizl99tSVJK1HnTP9S4FhmPpCZjwK3AVe0F8jMb2XmV4AfL9p2B3BnZp7MzIeAO4GdPai3JGkF6oT+BPBg2/JMta6O1WwrSeqxc2qUiQ7rsub+a20bEdcB1wFs2rSp5q4lqZ79h2fZN3WU43PzbBwbZc+ObezaXub5Z50z/Rng4rbli4DjNfdfa9vMvCUzJzNzcnx8vOauJWl5+w/PsveOe5mdmyeB2bl59t5xL/sPzw67akNRJ/QPAlsjYktEnAvsBg7U3P8UcFlEnF9dwL2sWidJA7Fv6ijzp848Zt38qTPsmzo6pBoN17Khn5mngetphfX9wKcy80hE3BgRlwNExIsiYgZ4PfDBiDhSbXsSeBetPxwHgRurdZI0EMfn5rta33R1xvTJzM8Bn1u07k/bHh+kNXTTadtbgVtXUUdJWrGNY6PMdgj4jWOjQ6jN8PmJXEmNtmfHNkY3jDxm3eiGEfbs2DakGg1XrTN9SVqvFmbpOHunxdCX1Hi7tk8UG/KLObwjSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVxNCXpIIY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFcTQl6SCGPqSVBBDX5IKYuhLUkEMfUkqiKEvSQUx9CWpIIa+JBXE0Jekghj6klQQQ1+SCmLoS1JBDH1JKoihL0kFMfQlqSCGviQVpFboR8TOiDgaEcci4oYOz58XEZ+snr8nIjZX6zdHxHxEfLn6+ZveVl+S1I1zlisQESPAzcCrgRngYEQcyMz72opdCzyUmb8QEbuB9wBXVc99MzOf3+N6S5JWoM6Z/qXAscx8IDMfBW4DrlhU5grgI9Xj24FXRkT0rpqSpF6oE/oTwINtyzPVuo5lMvM08DDw9Oq5LRFxOCL+LSJ+tdMBIuK6iJiOiOkTJ0501QBJUn11Qr/TGXvWLPNdYFNmbgfeDnw8Ip76uIKZt2TmZGZOjo+P16iSJGkl6oT+DHBx2/JFwPGlykTEOcDTgJOZ+aPM/CFAZh4Cvgk8a7WVliStTJ3QPwhsjYgtEXEusBs4sKjMAeCa6vGVwF2ZmRExXl0IJiKeCWwFHuhN1SVJ3Vp29k5mno6I64EpYAS4NTOPRMSNwHRmHgA+DHwsIo4BJ2n9YQB4OXBjRJwGzgBvysyT/WiIJGl5kbl4eH64Jicnc3p6etjVkKR1JSIOZebkcuX8RK4kFcTQl6SCLDumr/Vt/+FZ9k0d5fjcPBvHRtmzYxu7ti/+mMX60sQ2SYNi6DfY/sOz7L3jXuZPnQFgdm6evXfcC7BuQ7KJbZIGyeGdBts3dfQn4bhg/tQZ9k0dHVKNVq+JbZIGydBvsONz812tXw+a2CZpkAz9Bts4NtrV+vWgiW2SBsnQb7A9O7YxumHkMetGN4ywZ8e2IdVo9ZrYJmmQvJDbYAsXNps006WJbZIGyU/kSlID+IlcSdLjGPqSVBBDX5IK4oVcFcVbOKh0hr6K4S0cJENfPTCIs+deHONst3Aw9FUKQ1+rMoiz514dw1s4SF7I1SoN4gZovTqGt3CQDH2t0iDOnnt1jCbewmH/4Vle9u672HLDP/Kyd9/F/sOzw66S1jhDX6syiLPnXh1j1/YJbnrdc5kYGyWAibFRbnrdc9fteP7CsNfs3DzJT4e9DH6djWP6WpU9O7Y9Zrwden/23Mtj7No+sW5DfjEvTGslDH2tyiBugNaUm6z1epaTF6a1Eoa+Vm0QZ8/r/Qy9H7OcNo6NMtsh4L0wrbNxTF8agH7McmrihWn1n2f60gD0YyimKcNeGixD/yy8T4t6pV9DMet92EuD5/DOEpwOp15yKEZrhaG/hEF80lTlaNpnBLR+ObyzBKfDqdeaNhTj8Of65Jn+ErxPi7Q0hz/XL0N/CY7Bar0Z5H14ejn86f2DBsvhnSU4HU7ryaC/IKZXw59+sc3gGfpn0bQxWDXXoO/D06spqN4/aPAc3pEaYNATD3o1/NnUCRNrecjK0JcaYNATD3o1BbWJEybW+kVuh3fUU3Wm8fVrql/JUwgHcYvrxXox/NmLeq+1fl/rQ1ZFhf5ae3M0TZ2Lcv26cFf6BcH1OvFgtfVei/2+1oesaoV+ROwE3g+MAB/KzHcvev484KPAC4EfAldl5req5/YC1wJngLdk5lTPat+FtfjmaJo6Zzj9Ogta62dXg7BeJx6spt5rsd/X+i2vIzPPXiBiBPg68GpgBjgIXJ2Z97WV+QPgeZn5pojYDfx2Zl4VEZcAnwAuBTYC/wI8KzPPLD7OgsnJyZyenu66Ie1n8U8b3UAEzD1y6iePH3rk1JLbjnUoP/fIqcecdSy3/07b/sazx7n7aydqb9OPx0u1oZf1WzjG2z75ZTq9mwJ471XPZ9/U0Y6/DN30w1L9Xne/q+nT1fZ1neGu1fTvUvvs13uxV6/l4jZ0s6+6v9e9rOtyr+VDj5wi4DG/CwvLS73H677fzyYiDmXm5LLlaoT+S4E/y8wd1fJegMy8qa3MVFXmPyLiHOB7wDhwQ3vZ9nJLHW8lob/4LL6XRjeM8DsvnOAzh2b7sv9BGEQbRjeM8MQNT+j4Szg2uoEfnf7xqo49umHkcRcKV9rvw+zTxe3oxXu3H/vs5ti9eC0X2gA05nd5IegX/wGoo9P7fdnj1Qz9OrN3JoAH25ZnqnUdy2TmaeBh4Ok1t121Tv/F65X5U2f4xD0PrtvAh8G0Yf7UGTLpOI0vglUfu9OnPVfa78Ps08Xt6MV7tx/77ObYvXgtF9rQpN/lBEYiug586O/NHeuEfnRYt7gdS5Wpsy0RcV1ETEfE9IkTJ2pU6bH6fYHkzDL/G1oPBtGGh+dPdZzGN3eW/4J3Y3E/r6bfh9mn7fXu1Xu3H/usq1ev5fG5+cb9Lq/meP16LepcyJ0BLm5bvgg4vkSZmWp452nAyZrbkpm3ALdAa3inbuUXLHXhpFdGItZ98A+iDRvHRjtelFtqzL3bOi2+ELZUv9fZ7zD7tL0dvXrv9mOfdfXqtVxow0rrvhb7fTXH69eF3zpn+geBrRGxJSLOBXYDBxaVOQBcUz2+ErgrWxcLDgC7I+K8iNgCbAW+2Juq/1SnTwcuZXTDCG94yaauyl/94otrl1+LBtGGs82tXurTm93UqdP+V7rfYfbp4nZ0894d5D67OXYvXsuFNqy07mux31dzvH5+xmLZM/3MPB0R1wNTtKZs3pqZRyLiRmA6Mw8AHwY+FhHHaJ3h7662PRIRnwLuA04Dbz7bzJ2VWjzXt85V8clnXNDz8mt99k57G/oxe2epi05nm4vdbT+sdr+r6dN+zN6p+95dzT4HMXtnNa9lpzZ0+7s2qH5fSV+vJDv6+RmLZWfvDNpKp2xKUsl6OXtHktQQhr4kFcTQl6SCGPqSVBBDX5IKsuZm70TECeDbq9jFhcAPelSd9aLENkOZ7S6xzVBmu7tt8zMyc3y5Qmsu9FcrIqbrTFtqkhLbDGW2u8Q2Q5nt7lebHd6RpIIY+pJUkCaG/i3DrsAQlNhmKLPdJbYZymx3X9rcuDF9SdLSmnimL0laQmNCPyJ2RsTRiDgWETcMuz79EhEXR8TdEXF/RByJiLdW6y+IiDsj4hvVv+cPu669FhEjEXE4Ij5bLW+JiHuqNn+yuvV3o0TEWETcHhFfq/r8pU3v64h4W/Xe/mpEfCIintjEvo6IWyPi+xHx1bZ1Hfs2Wv66yrevRMQLVnrcRoR+9eXtNwOvAS4Brq6+lL2JTgN/lJnPAV4CvLlq6w3A5zNzK/D5arlp3grc37b8HuC9VZsfAq4dSq366/3AP2fms4FfotX+xvZ1REwAbwEmM/MXad3OfTfN7Ou/A3YuWrdU376G1veRbAWuAz6w0oM2IvSBS4FjmflAZj4K3AZcMeQ69UVmfjczv1Q9/l9aITBBq70fqYp9BNg1nBr2R0RcBPwm8KFqOYBXALdXRZrY5qcCL6f1fRVk5qOZOUfD+5rW93yMVt/C9yTguzSwrzPz32l9/0i7pfr2CuCj2fIFYCwifm4lx21K6A/kC9jXmojYDGwH7gF+NjO/C60/DMDPDK9mffE+4E+AH1fLTwfmMvN0tdzEPn8mcAL422pY60MR8WQa3NeZOQv8JfAdWmH/MHCI5vf1gqX6tmcZ15TQr/UF7E0SEU8BPgP8YWb+z7Dr008R8VvA9zPzUPvqDkWb1ufnAC8APpCZ24H/o0FDOZ1UY9hXAFuAjcCTaQ1tLNa0vl5Oz97vTQn9Wl/A3hQRsYFW4P99Zt5Rrf7vhf/uVf9+f1j164OXAZdHxLdoDd29gtaZ/1g1BADN7PMZYCYz76mWb6f1R6DJff0q4L8y80RmngLuAH6Z5vf1gqX6tmcZ15TQr/Pl7Y1QjWV/GLg/M/+q7an2L6e/BviHQdetXzJzb2ZelJmbafXtXZn5e8DdwJVVsUa1GSAzvwc8GBEL35D9SlrfN93YvqY1rPOSiHhS9V5faHOj+7rNUn17APj9ahbPS4CHF4aBupaZjfgBXgt8Hfgm8I5h16eP7fwVWv+t+wrw5erntbTGuD8PfKP694Jh17VP7f914LPV42cCXwSOAZ8Gzht2/frQ3ucD01V/7wfOb3pfA38OfA34KvAx4Lwm9jXwCVrXLU7ROpO/dqm+pTW8c3OVb/fSmt20ouP6iVxJKkhThnckSTUY+pJUEENfkgpi6EtSQQx9SSqIoS9JBTH0Jakghr4kFeT/AXXO95k7YHX+AAAAAElFTkSuQmCC\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stdout",
"output_type": "stream",
"text": [
" HMCR: status=ProblemStatus.PrimalAndDualFeasible time=0.64s\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAX0AAAD8CAYAAACb4nSYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAE5xJREFUeJzt3W+MXNd93vHv0xUlr2PElK1tUZKSSSEMY7VuzXQiO3XrtoltymkhEakCy6hRBRAgtIjQtG5ZiPALt0oBx2HRf4CQSojdukFj2VYIlgiQEKqktG8qhUvTFU3JrCnVkUg6FVOZblEvLJL+9cVe2sP1rnaGnOXszPl+gAXnnntm5nf2DJ+dOXNnbqoKSVIb/sS4C5AkXTuGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakh1427gKVuuumm2rp167jLkKSJcuTIkT+uqrnV+q270N+6dSvz8/PjLkOSJkqSPxykn8s7ktQQQ1+SGmLoS1JDDH1JashAoZ/kjiQnkpxM8uAy+z+e5PkkzyV5Msk7+vZdTPKV7ufgKIuXJA1n1aN3kswADwMfBE4Bh5McrKrn+7odBXpV9Z0kfxf4NeAj3b6Fqnr3iOuWJtqBo6fZd+gEZ84tsGnjLHt27WD3zs3jLksNGOSZ/u3Ayap6qapeBx4D7urvUFVPV9V3us1ngC2jLVOaHgeOnmbv/mOcPrdAAafPLbB3/zEOHD097tLUgEFCfzPwSt/2qa5tJfcBv9u3/aYk80meSbL7CmqUpsq+QydYOH/xsraF8xfZd+jEmCpSSwb5cFaWaVv2xLpJPgb0gL/S13xLVZ1JcivwVJJjVfXikuvdD9wPcMsttwxUuDSpzpxbGKpdGqVBnumfAm7u294CnFnaKckHgE8Ad1bVdy+1V9WZ7t+XgN8Hdi69blU9WlW9qurNza36KWJpom3aODtUuzRKg4T+YWB7km1JrgfuAS47CifJTuARFgP/1b72G5Pc0F2+CXgf0P8GsNScPbt2MLth5rK22Q0z7Nm1Y0wVqSWrLu9U1YUkDwCHgBngs1V1PMlDwHxVHQT2AW8BvpQE4OWquhN4J/BIku+x+AfmV5cc9SM159JROh69o3FI1bLL82PT6/XKL1yTpOEkOVJVvdX6+YlcSWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaMlDoJ7kjyYkkJ5M8uMz+jyd5PslzSZ5M8o6+ffcm+Xr3c+8oi5ckDWfV0E8yAzwMfBi4DfhoktuWdDsK9KrqzwGPA7/WXfdtwCeB9wC3A59McuPoypckDWOQZ/q3Ayer6qWqeh14DLirv0NVPV1V3+k2nwG2dJd3AU9U1WtV9S3gCeCO0ZQuSRrWIKG/GXilb/tU17aS+4DfHea6Se5PMp9k/uzZswOUJEm6EoOEfpZpq2U7Jh8DesC+Ya5bVY9WVa+qenNzcwOUJEm6EoOE/ing5r7tLcCZpZ2SfAD4BHBnVX13mOtKkq6NQUL/MLA9ybYk1wP3AAf7OyTZCTzCYuC/2rfrEPChJDd2b+B+qGuTJI3Bdat1qKoLSR5gMaxngM9W1fEkDwHzVXWQxeWctwBfSgLwclXdWVWvJfkVFv9wADxUVa+tyUgkSatK1bLL82PT6/Vqfn5+3GVI0kRJcqSqeqv18xO5ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIas+olcSevXgaOn2XfoBGfOLbBp4yx7du1g9843+hJctc7QlybUgaOn2bv/GAvnLwJw+twCe/cfAzD4tSKXd6QJte/Qie8H/iUL5y+y79CJMVWkSWDoSxPqzLmFodolMPSlibVp4+xQ7RIY+tLE2rNrB7MbZi5rm90ww55dO8ZUkSaBb+RKE+rSm7UevaNhGPrSBNu9c7Mhr6G4vCNJDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGjJQ6Ce5I8mJJCeTPLjM/vcn+XKSC0nuXrLvYpKvdD8HR1W4JGl4q36ffpIZ4GHgg8Ap4HCSg1X1fF+3l4FfBP7RMjexUFXvHkGtkqSrNMhJVG4HTlbVSwBJHgPuAr4f+lX1jW7f99agRknSiAyyvLMZeKVv+1TXNqg3JZlP8kyS3UNVJ0kaqUGe6WeZthriPm6pqjNJbgWeSnKsql687A6S+4H7AW655ZYhblqSNIxBnumfAm7u294CnBn0DqrqTPfvS8DvAzuX6fNoVfWqqjc3NzfoTUuShjRI6B8GtifZluR64B5goKNwktyY5Ibu8k3A++h7L0CSdG2tGvpVdQF4ADgEvAB8saqOJ3koyZ0ASX4qySngF4BHkhzvrv5OYD7JfweeBn51yVE/kqRrKFXDLM+vvV6vV/Pz8+MuQ5ImSpIjVdVbrZ+fyJWkhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGDfPeOtK4cOHqafYdOcObcAps2zrJn1w527xzmOwCldhn6migHjp5m7/5jLJy/CMDpcwvs3X8MwOCXBuDyjibKvkMnvh/4lyycv8i+QyfGVJE0WQx9TZQz5xaGapd0OUNfE2XTxtmh2iVdztDXRNmzawezG2Yua5vdMMOeXTvGVJE0WXwjVxPl0pu1Hr0jXRlDXxNn987Nhrx0hVzekaSGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1ZKDQT3JHkhNJTiZ5cJn970/y5SQXkty9ZN+9Sb7e/dw7qsIlScNbNfSTzAAPAx8GbgM+muS2Jd1eBn4R+K0l130b8EngPcDtwCeT3Hj1ZUuSrsQgz/RvB05W1UtV9TrwGHBXf4eq+kZVPQd8b8l1dwFPVNVrVfUt4AngjhHULUm6AoOE/mbglb7tU13bIK7mupKkERsk9LNMWw14+wNdN8n9SeaTzJ89e3bAm5YkDWuQ0D8F3Ny3vQU4M+DtD3Tdqnq0qnpV1ZubmxvwpiVJwxok9A8D25NsS3I9cA9wcMDbPwR8KMmN3Ru4H+raJEljsGroV9UF4AEWw/oF4ItVdTzJQ0nuBEjyU0lOAb8APJLkeHfd14BfYfEPx2Hgoa5NkjQGqRp0ef7a6PV6NT8/P+4yJGmiJDlSVb3V+vmJXElqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDDH1JaoihL0kNMfQlqSGGviQ1xNCXpIYY+pLUEENfkhpi6EtSQwx9SWqIoS9JDTH0Jakhhr4kNcTQl6SGGPqS1BBDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGjJQ6Ce5I8mJJCeTPLjM/huSfKHb/2ySrV371iQLSb7S/fzb0ZYvSRrGdat1SDIDPAx8EDgFHE5ysKqe7+t2H/CtqvqxJPcAnwY+0u17sarePeK6JUlXYJBn+rcDJ6vqpap6HXgMuGtJn7uAz3WXHwd+NklGV6YkaRQGCf3NwCt926e6tmX7VNUF4NvA27t925IcTfJfkvzl5e4gyf1J5pPMnz17dqgBSJIGN0joL/eMvQbs803glqraCXwc+K0kP/pDHaserapeVfXm5uYGKEmSdCUGCf1TwM1921uAMyv1SXId8Fbgtar6blX9b4CqOgK8CPz41RYtSboyg4T+YWB7km1JrgfuAQ4u6XMQuLe7fDfwVFVVkrnujWCS3ApsB14aTemSpGGtevROVV1I8gBwCJgBPltVx5M8BMxX1UHgM8BvJjkJvMbiHwaA9wMPJbkAXAT+TlW9thYDkSStLlVLl+fHq9fr1fz8/LjLkKSJkuRIVfVW6+cnciWpIasu70gajQNHT7Pv0AnOnFtg08ZZ9uzawe6dS49+ltaWoS9dAweOnmbv/mMsnL8IwOlzC+zdfwzA4Nc15fKOdA3sO3Ti+4F/ycL5i+w7dGJMFalVPtPXG3JJYjTOnFsYql1aKz7T14ouLUmcPrdA8YMliQNHT4+7tImzaePsUO3SWjH0tSKXJEZnz64dzG6YuaxtdsMMe3btGFNFapXLO1qRSxKjc2lJzKUyjZuhrxVt2jjL6WUC3iWJK7N752ZDXmPn8o5W5JKENH2aeqbvkSjDcUlCmj7NhL4fjrkyLklI06WZ0H+jI1EMtfXPV2nSaDQT+h6JMrl8lSaNTjNv5PrhmMnl5wWk0Wkm9D0SZXL5Kk0anWZCf/fOzXzq59/F5o2zBNi8cZZP/fy7XB6YAL5Kk0anmTV98EiUSbVn147L1vTBV2nSlWoq9DWZ/LyANDqGviaCr9Kk0WhmTV+SZOhLUlMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhhj6ktQQQ1+SGmLoS1JDBgr9JHckOZHkZJIHl9l/Q5IvdPufTbK1b9/erv1Ekl2jK12SNKxVv3snyQzwMPBB4BRwOMnBqnq+r9t9wLeq6seS3AN8GvhIktuAe4A/A2wC/nOSH6+qy8+IMQL9p9N76+wGEjj3nfOXXe7/oq5B+l/N5U0bZ/lrPzHH0187u2b3cS3qu5a/s0FquJJ5X4+/y1GNYZD5We+PxUmq9VrUt9ZfKJiqeuMOyU8D/6SqdnXbewGq6lN9fQ51ff5bkuuAPwLmgAf7+/b3W+n+er1ezc/PDzWIpafTeyOzG2b4m39hM7995PRA/bU+fmezG2Z+6PwHw8z7erF0HKMYw3qYH43Wco/31SQ5UlW91foNsryzGXilb/tU17Zsn6q6AHwbePuA171qy51ObyUL5y/y+Wdf8T/HENbD72y50yMOM+/rxdJxjGIM62F+NFpreTrQQUI/y7QtfXmwUp9BrkuS+5PMJ5k/e/bsACVdbtjT5l1c5dWNfth6+J0tnedJPV1if92jGsN6mB+N1lo9vgcJ/VPAzX3bW4AzK/XplnfeCrw24HWpqkerqldVvbm5ucGr7wx72ryZLPe3SG9kPfzOls7zpJ4usb/uUY1hPcyPRmutHt+DhP5hYHuSbUmuZ/GN2YNL+hwE7u0u3w08VYtvFhwE7umO7tkGbAf+YDSl/8ByJz1fyeyGGT76npsH7q/18Ttb7vSIw8z7erF0HKMYw3qYH43WWp4OdNXQ79boHwAOAS8AX6yq40keSnJn1+0zwNuTnAQ+zg/ewD0OfBF4Hvg94JfW4sidpSc93zi7gRvfvOGHLl86Gfo/2/2ugfpfzeXNG2f52HtvWdP7uBb1Xcvf2Wo1LH1Ta9B5X2+/y/5xXO0YBpmf9f5YnKRar0V9Kz3eR2XVo3eutSs5ekeSWjfKo3ckSVPC0Jekhhj6ktQQQ1+SGmLoS1JD1t3RO0nOAn94FTdxE/DHIypnUrQ4Zmhz3C2OGdoc97BjfkdVrfrp1nUX+lcryfwghy1NkxbHDG2Ou8UxQ5vjXqsxu7wjSQ0x9CWpIdMY+o+Ou4AxaHHM0Oa4WxwztDnuNRnz1K3pS5JWNo3P9CVJK5ia0F/t5O3TIsnNSZ5O8kKS40l+uWt/W5Inkny9+/fGcdc6aklmkhxN8jvd9rYkz3Zj/kL31d9TJcnGJI8n+Vo35z897XOd5B90j+2vJvl8kjdN41wn+WySV5N8ta9t2bnNon/T5dtzSX7ySu93KkK/7+TtHwZuAz7anZR9Gl0A/mFVvRN4L/BL3VgfBJ6squ3Ak932tPllFr/e+5JPA/+yG/O3gPvGUtXa+tfA71XVTwB/nsXxT+1cJ9kM/D2gV1V/Fphh8Rwe0zjX/x64Y0nbSnP7YRbPR7IduB/49Su906kIfeB24GRVvVRVrwOPAXeNuaY1UVXfrKovd5f/L4shsJnF8X6u6/Y5YPd4KlwbSbYAfx34jW47wM8Aj3ddpnHMPwq8n8XzVVBVr1fVOaZ8roHrgNnuLHxvBr7JFM51Vf1XFs8w2G+lub0L+A+16BlgY5I/fSX3Oy2hf01OwL7eJNkK7ASeBf5UVX0TFv8wAH9yfJWtiX8F/GPge93224Fz3Ul+YDrn/FbgLPDvumWt30jyI0zxXFfVaeCfAy+zGPbfBo4w/XN9yUpzO7KMm5bQH+gE7NMkyVuA3wb+flX9n3HXs5aS/A3g1ao60t+8TNdpm/PrgJ8Efr2qdgL/jylayllOt4Z9F7AN2AT8CItLG0tN21yvZmSP92kJ/YFOwD4tkmxgMfD/Y1Xt75r/16WXe92/r46rvjXwPuDOJN9gcenuZ1h85r+xWwKA6ZzzU8Cpqnq2236cxT8C0zzXHwD+Z1WdrarzwH7gLzL9c33JSnM7soybltAf5OTtU6Fby/4M8EJV/Yu+Xf0np78X+E/Xura1UlV7q2pLVW1lcW6fqqq/BTwN3N11m6oxA1TVHwGvJLl0huyfZfF801M71ywu67w3yZu7x/qlMU/1XPdZaW4PAn+7O4rnvcC3Ly0DDa2qpuIH+DngfwAvAp8Ydz1rOM6/xOLLuueAr3Q/P8fiGveTwNe7f9827lrXaPx/Ffid7vKtwB8AJ4EvATeMu741GO+7gfluvg8AN077XAP/FPga8FXgN4EbpnGugc+z+L7FeRafyd+30tyyuLzzcJdvx1g8uumK7tdP5EpSQ6ZleUeSNABDX5IaYuhLUkMMfUlqiKEvSQ0x9CWpIYa+JDXE0Jekhvx/f7KT32zcZlYAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"n, m = 100,2000\n",
"alpha = 0.9\n",
"rbar = 0.005\n",
"p = [1.0/m]*m # Uniform probability\n",
"\n",
"MCVaR = modelCVaR(alpha, rbar, p, r[0:n,0:m], m, n)\n",
"MLogExpCR = modelLogExpCR(alpha, rbar, p, r[0:n,0:m], m, n)\n",
"MHMCR = modelHMCR(alpha, rbar, p, r[0:n,0:m], m, n, q=5.9)\n",
"\n",
"for M in [MCVaR, MLogExpCR, MHMCR]:\n",
" M.setSolverParam(\"numThreads\", 4)\n",
" M.solve()\n",
" print('{0: >8}: status={1} time={2:.2f}s'.format(M.getName(), M.getProblemStatus(), M.getSolverDoubleInfo(\"optimizerTime\")))\n",
" plt.plot(M.getVariable('X').level(), 'o')\n",
" plt.show()\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"
This work is licensed under a Creative Commons Attribution 4.0 International License. The **MOSEK** logo and name are trademarks of Mosek ApS. The code is provided as-is. Compatibility with future release of **MOSEK** or the `Fusion API` are not guaranteed. For more information contact our [support](mailto:support@mosek.com). "
]
}
],
"metadata": {
"anaconda-cloud": {},
"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.6.4"
}
},
"nbformat": 4,
"nbformat_minor": 1
}