{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "![MOSEK ApS](https://www.mosek.com/static/images/branding/webgraphmoseklogocolor.png )" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Optimization problem over nonnegative trigonometric polynomials. #\n", "\n", "We consider the trigonometric polynomial\n", " \n", "$$H(\\omega) = x_0 + 2 \\sum_{k=1}^n [ \\Re(x_k) \\cos(\\omega k) + \\Im(x_k) \\sin(\\omega k) ],$$\n", "\n", "where $H(\\omega)$ is a real valued function paramtrized by the complex vector $x\\in {\\bf C}^{n+1}$, and where $\\Im(x_0)=0$.\n", "\n", "The example shows how to construct a *non-negative* polynomial $H(\\omega)\\geq 0, \\: \\forall \\omega$ that satisfies,\n", " \n", "$$1 - \\delta \\leq H(\\omega) \\leq 1 + \\delta, \\quad \\forall \\omega \\in [0, \\omega_p]$$\n", "\n", "while minimizing $\\sup_{\\omega\\in [\\omega_s,\\pi]} H(\\omega)$ over the variables $x$.\n", "\n", "In the signal processing literature such a trigonometric polynomial is known as (the squared amplitude response of) a Chebyshev lowpass filter. \n", "\n", "A squared amplitude response $H(\\omega)$ is always symmetric around $0$, so $\\Im(x_k)=0$, and we consider only\n", "\n", "$$H(\\omega) = x_0 + 2 \\sum_{k=1}^n x_k \\cos(\\omega k) $$\n", "\n", "over the real vector $x\\in {\\bf R}^{n+1}$. However, the concepts in the example are readily applied to the case with $x\\in {\\bf C}^{n+1}$.\n", "\n", "References:\n", " 1. \"Squared Functional Systems and Optimization Problems\", Y. Nesterov, 2004.\n", " \n", " 2. \"Convex Optimization of Non-negative Polynomials: Structured Algorithms and Applications\", Ph.D thesis, Y. Hachez, 2003.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import mosek\n", "from mosek.fusion import *\n", "from math import cos, pi\n", "import numpy as np\n", "import sys" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Nonnegativity everywhere ###\n", "\n", "Nesterov proved in [1] that $H(\\omega) \\geq 0, \\: \\forall \\omega$ if and only if \n", "$$x_i = \\langle T_i^{n+1}, X \\rangle, \\quad X \\in {\\mathcal H}^{n+1}_+,$$\n", "where ${\\mathcal H}_+$ is the cone of Hermitian semidefinite matrices and $T_i$ is a Toeplitz matrix\n", "$$[T_i]_{kl} = \\left \\{ \\begin{array}{ll}1, & k-l=i\\\\0 & \\text{otherwise}.\\end{array} \\right .$$\n", "For example, for $n=2$ we have\n", "$$\n", " T_0 = \\left[\\begin{array}{ccc}\n", " 1 & 0 & 0\\\\0 & 1 & 0\\\\0 & 0 & 1\n", " \\end{array}\n", " \\right], \\quad\n", " T_1 = \\left[\\begin{array}{ccc}\n", " 0 & 0 & 0\\\\1 & 0 & 0\\\\0 & 1 & 0\n", " \\end{array}\n", " \\right],\n", " \\quad\n", " T_2 = \\left[\\begin{array}{ccc}\n", " 0 & 0 & 0\\\\0 & 0 & 0\\\\1 & 0 & 0\n", " \\end{array}\n", " \\right].\n", "$$\n", "In our case we have $\\Im(x_i)=0$, i.e., $X\\in {\\mathcal S}^{n+1}_+$ is a real symmetric semidefinite matrix.\n", "\n", "We define the *cone on nonnegative trigonometric polynomials* as\n", "$$\n", " K^n_{0,\\pi} = \\{ x\\in \\mathbf{R} \\times \\mathbf{C}^n \\mid x_i = \\langle X, T_i\\rangle, \\: i=0,\\dots,n, \\: X\\in\\mathcal{H}_+^{n+1} \\}.\n", "$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def T_dot_X(n, i, X, a=1.0):\n", " if i>=n or i<=-n: return Expr.constTerm(0.0)\n", " else: return Expr.dot(Matrix.diag(n, a, -i), X)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def trigpoly_0_pi(M, x):\n", " '''x[i] == '''\n", " n = int(x.getSize()-1)\n", " X = M.variable(\"X\", Domain.inPSDCone(n+1))\n", " \n", " for i in range(n+1):\n", " M.constraint(Expr.sub(T_dot_X(n+1,i,X), x.index(i)), Domain.equalsTo(0.0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we have dropped the imaginary part of $X$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Nonnegativity on $[0,a]$ ###\n", "\n", "Similarly, $H(\\omega)$ is nonnegative on $[0,a]$ if and only if\n", "\n", "$$x_i =\n", "\\langle X_1, T_i^{n+1} \\rangle + \n", "\\langle X_2, T_{i+1}^n \\rangle +\n", "\\langle X_2, T_{i-1}^n \\rangle -\n", "2 \\cos(a)\\langle X_2, T_{i}^n \\rangle, \\quad \n", "X_1 \\in {\\mathcal H}^{n+1}_+, \\:\n", "X_2 \\in {\\mathcal H}^n_+,\n", "$$\n", "or equivalently\n", "$$\n", " K^n_{0,a} = \\{ x\\in \\mathbf{R}\\times \\mathbf{C}^n \\mid\n", " x_i = \\langle X_1, T_i^{n+1} \\rangle +\n", " \\langle X_2 , T_{i+1}^n \\rangle +\n", " \\langle X_2 , T_{i-1}^n \\rangle -\n", " 2\\cos(a)\\langle X_2 , T_i^n\\rangle, \\: X_1\\in \\mathcal{H}_+^{n+1}, X_2\\in\\mathcal{H}_+^n \\}.\n", "$$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def trigpoly_0_a(M, x, a):\n", " '''x[i] == + + - 2*cos(a)*'''\n", " n = int(x.getSize()-1)\n", " X1 = M.variable(Domain.inPSDCone(n+1))\n", " X2 = M.variable(Domain.inPSDCone(n))\n", "\n", " for i in range(n+1):\n", " M.constraint(Expr.sub(Expr.add([ T_dot_X(n+1,i,X1), \n", " T_dot_X(n,i+1,X2), \n", " T_dot_X(n,i-1,X2), \n", " T_dot_X(n,i,X2,-2*cos(a))]),\n", " x.index(i)), Domain.equalsTo(0.0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we have dropped the imaginary part of $X_1$ and $X_2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Nonnegativity on $[a,\\pi]$ ###\n", "\n", "Finally, $H(\\omega)$ is nonnegative on $[a,\\pi]$ if and only if\n", "\n", "$$x_i = \n", "\\langle X_1, T_i^{n+1} \\rangle -\n", "\\langle X_2, T_{i+1}^n \\rangle -\n", "\\langle X_2, T_{i-1}^n \\rangle +\n", "2 \\cos(a)\\langle X_2, T_{i}^n \\rangle, \\quad \n", "X_1 \\in {\\mathcal S}^{n+1}_+, \\:\n", "X_2 \\in {\\mathcal S}^n_+,\n", "$$\n", "or equivalently\n", "$$\n", " K^n_{a,\\pi} = \\{ x\\in \\mathbf{R}\\times \\mathbf{C}^n \\mid\n", " x_i = \\langle X_1, T_i^{n+1} \\rangle -\n", " \\langle X_2 , T_{i+1}^n \\rangle -\n", " \\langle X_2 , T_{i-1}^n \\rangle +\n", " 2\\cos(a)\\langle X_2 , T_i^n\\rangle, \\: X_1\\in \\mathcal{H}_+^{n+1}, X_2\\in\\mathcal{H}_+^n \\}.\n", "$$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def trigpoly_a_pi(M, x, a):\n", " '''x[i] == - - + 2*cos(a)*'''\n", " n = int(x.getSize()-1)\n", " X1 = M.variable(Domain.inPSDCone(n+1))\n", " X2 = M.variable(Domain.inPSDCone(n))\n", "\n", " for i in range(n+1):\n", " M.constraint(Expr.add([ T_dot_X(n+1,i,X1), \n", " T_dot_X(n,i+1,X2,-1.0), \n", " T_dot_X(n,i-1,X2,-1.0), \n", " T_dot_X(n,i,X2,2*cos(a)),\n", " Expr.neg(x.index(i)) ]),\n", " Domain.equalsTo(0.0))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we have dropped the imaginary part of $X_1$ and $X_2$." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## An epigraph formulation ##\n", "The *epigraph* $H(\\omega) \\leq t$ can now be characterized simply as \n", "$-u \\in K^n_{[a,b]}, \\: u=(x_0-t, \\, x_{1:n}).$" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "def epigraph(M, x, t, a, b):\n", " '''Models 0 <= H(w) <= t, for all w in [a, b], where\n", " H(w) = x0 + 2*x1*cos(w) + 2*x2*cos(2*w) + ... + 2*xn*cos(n*w)'''\n", " n = int(x.getSize()-1) \n", " u = M.variable(n+1, Domain.unbounded())\n", " M.constraint(Expr.sub(t,Expr.add(x.index(0), u.index(0))), Domain.equalsTo(0.0)) \n", " M.constraint(Expr.add(x.slice(1,n+1), u.slice(1,n+1)), Domain.equalsTo(0.0))\n", " \n", " if a==0.0 and b==pi:\n", " trigpoly_0_pi(M, u)\n", " elif a==0.0 and b= 0\n", "trigpoly_0_pi(M, x)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Passband specifications ###" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "wp = pi/4\n", "delta = 0.05\n", "\n", "# H(w) <= (1+delta), w in [0, wp] \n", "epigraph(M, x, 1.0+delta, 0.0, wp)\n", "\n", "# (1-delta) <= H(w), w in [0, wp]\n", "hypograph(M, x, 1.0-delta, 0.0, wp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Stopband specifications ###" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "ws = wp + pi/8\n", "\n", "# H(w) < t, w in [ws, pi]\n", "t = M.variable(\"t\", 1, Domain.greaterThan(0.0))\n", "epigraph(M, x, t, ws, pi)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Setting the objective ###" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "M.objective(ObjectiveSense.Minimize, t)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "M.setLogHandler(sys.stdout)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Problem\n", " Name : trigpoly \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 67 \n", " Cones : 0 \n", " Scalar variables : 35 \n", " Matrix variables : 7 \n", " Integer variables : 0 \n", "\n", "Optimizer started.\n", "Presolve started.\n", "Linear dependency checker started.\n", "Linear dependency checker terminated.\n", "Eliminator started.\n", "Freed constraints in eliminator : 23\n", "Eliminator terminated.\n", "Eliminator started.\n", "Freed constraints in eliminator : 0\n", "Eliminator terminated.\n", "Eliminator - tries : 2 time : 0.00 \n", "Lin. dep. - tries : 1 time : 0.00 \n", "Lin. dep. - number : 0 \n", "Presolve terminated. Time: 0.01 \n", "Problem\n", " Name : trigpoly \n", " Objective sense : min \n", " Type : CONIC (conic optimization problem)\n", " Constraints : 67 \n", " Cones : 0 \n", " Scalar variables : 35 \n", " Matrix variables : 7 \n", " Integer variables : 0 \n", "\n", "Optimizer - threads : 20 \n", "Optimizer - solved problem : the primal \n", "Optimizer - Constraints : 44\n", "Optimizer - Cones : 1\n", "Optimizer - Scalar variables : 13 conic : 12 \n", "Optimizer - Semi-definite variables: 7 scalarized : 429 \n", "Factor - setup time : 0.00 dense det. time : 0.00 \n", "Factor - ML order time : 0.00 GP order time : 0.00 \n", "Factor - nonzeros before factor : 990 after factor : 990 \n", "Factor - dense dim. : 0 flops : 1.18e+05 \n", "ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME \n", "0 1.9e+01 1.0e+00 1.0e+00 0.00e+00 0.000000000e+00 0.000000000e+00 1.0e+00 0.03 \n", "1 6.0e+00 3.2e-01 3.2e-01 1.56e+00 -6.168509433e-03 -1.704940542e-02 3.2e-01 0.04 \n", "2 1.5e+00 8.0e-02 8.0e-02 1.59e+00 -3.629686898e-02 -5.729238420e-02 8.0e-02 0.05 \n", "3 3.1e-01 1.7e-02 1.7e-02 1.64e+00 -4.704273187e-03 -7.156318767e-03 1.7e-02 0.05 \n", "4 1.0e-01 5.6e-03 5.6e-03 1.91e-01 1.247879338e-03 3.220303293e-03 5.6e-03 0.05 \n", "5 5.8e-02 3.1e-03 3.1e-03 4.28e-01 4.186115053e-02 3.996134158e-02 3.1e-03 0.05 \n", "6 1.1e-02 5.9e-04 5.9e-04 4.33e-01 6.583025991e-02 6.538915159e-02 5.9e-04 0.05 \n", "7 2.7e-03 1.4e-04 1.4e-04 9.75e-01 7.102893174e-02 7.092631673e-02 1.4e-04 0.05 \n", "8 4.6e-04 2.4e-05 2.4e-05 1.01e+00 7.246509927e-02 7.245345168e-02 2.4e-05 0.06 \n", "9 3.3e-05 1.8e-06 1.8e-06 1.00e+00 7.271267688e-02 7.271184530e-02 1.8e-06 0.06 \n", "10 3.6e-06 1.9e-07 1.9e-07 1.00e+00 7.272951727e-02 7.272942593e-02 1.9e-07 0.06 \n", "11 5.8e-07 3.1e-08 3.1e-08 1.00e+00 7.273126748e-02 7.273125303e-02 3.1e-08 0.06 \n", "12 7.0e-08 3.9e-09 3.7e-09 1.00e+00 7.273156244e-02 7.273156074e-02 3.7e-09 0.06 \n", "13 6.3e-09 1.9e-09 3.4e-10 1.00e+00 7.273159968e-02 7.273159953e-02 3.4e-10 0.07 \n", "Optimizer terminated. Time: 0.07 \n", "\n", "\n", "Interior-point solution summary\n", " Problem status : PRIMAL_AND_DUAL_FEASIBLE\n", " Solution status : OPTIMAL\n", " Primal. obj: 7.2731599679e-02 nrm: 1e+00 Viol. con: 6e-09 var: 0e+00 barvar: 0e+00 \n", " Dual. obj: 7.2731599528e-02 nrm: 2e+00 Viol. con: 0e+00 var: 3e-10 barvar: 2e-09 \n" ] } ], "source": [ "M.solve()" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 0.33580937, 0.25254317, 0.1388759 , 0.02044169, -0.04797007,\n", " -0.05156956, -0.01627207, 0.01837092, 0.02948826, 0.02271378,\n", " -0.00952669])" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "x.level()" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([0.0727316])" ] }, "execution_count": 16, "metadata": {}, "output_type": "execute_result" } ], "source": [ "t.level()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the amplitude response ###" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAD8CAYAAACMwORRAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzt3XtUVXXeP/D3h4viBVETBW+hqKnk5cBJpAywNPGSZk1Ol3Gq1TxUs2wyny5TM9rU5DyreSpb489Kp8dsctSc+lWOqWC/sryiIHhPRfOClyDwgooinM/vD5AHEWQD+5y9zznv11qsdS5f9n7vTr7P5rvPPltUFURE5FsCrA5ARETmY7kTEfkgljsRkQ9iuRMR+SCWOxGRD2K5ExH5IJY7EZEPYrkTEfkgljsRkQ8KsmrFHTp00KioKKtWT0TklbKysn5W1fD6xllW7lFRUcjMzLRq9UREXklEDhsZx2kZIiIfxHInIvJBLHciIh/Ecici8kEsdyIiH8RyJyLyQSx3IiIfZNnn3Btr/fr1+Oabb3DjjTciKioKvXv3RmRkpNWx3Or48eNYsWIFAgMDER0djejoaHTu3BkiYnU0IrIpryv3DRs2YMaMGVc9dtttt2Hy5MmYNGkS2rVrZ1Eyc5WWluL999/H4sWLsWnTpmueT0xMxEcffQSe5UtEtfG6aZnnn38eJSUl2Lt3L9LS0jBz5kwUFRXhySefRGRkJF5++WUUFxdbHbNJcnJyMGTIEDzzzDO4dOkSXn/9dezcuRP79+9HWloa3njjDWRnZ2PgwIH48MMPwYucE9E1VNWSn7i4ODWLy+XSrKwsffjhhxWARkRE6Pz587W8vNy0dXhCWVmZvvLKKxoUFKQRERH6xRdf1Dn2xx9/1KSkJAWgTz75pAdTEpGVAGSqgY71iXKvbuPGjRofH68AdMyYMVpQUOCW9ZjtwoULeu+99yoA/dWvfqWFhYX1/k55eblOmzZNAeiHH37o/pBEZDm/LXfVitKbPXu2NmvWTDt37qxr1qxx27rMUFhYqMOGDVMR0XfeeadBv1tWVqbDhw/XFi1a6I4dO9yUkIjswmi5e92cuxEBAQGYMmUKNm3ahFatWuGOO+7AW2+9Zcu56WPHjmHYsGHYvHkzlixZgmeeeaZBvx8YGIhFixahTZs2uP/++3Hu3Dk3JSUib1JvuYvIfBHJF5GddTwvIvI3EckVke0iEmt+zMZxOBzIysrCxIkT8dxzzyE1NRWlpaVWx6py7NgxDB8+HHl5eUhLS8OkSZMatZyIiAgsXrwY+/bta/CbAxH5JiN77gsApFzn+dEAelf+pAJ4r+mxzBMaGoqlS5fiD3/4Az744AOkpKSgqKjI6ljIy8tDcnIyTp48ibS0NCQnJzdpecOHD8ezzz6LBQsW4IcffjAnJBF5rXrLXVW/B3C9NpwA4B+V00GbALQVEVudVRQQEIDXX38d//jHP7B+/XoMGzYMhw8b+r57tzh69CiSk5Px008/IS0tDQkJCaYs98UXX0RISAhef/11U5ZHRN7LjDn3LgCOVrufV/mY7UyePBnp6ek4ceIEhg4diuzsbI9nOHLkCJKTk1FQUID09HTTih0AwsPDMWXKFCxevJh770R+zoxyr+0c+FqPXIpIqohkikhmQUGBCatuuKSkJKxbtw7BwcFITEzEqlWrPLbuw4cPIzk5GT///DPS09MxdOhQ09fx3HPPce+diEwp9zwA3ard7wrgeG0DVXWeqjpV1RkeXu/1Xd0mJiYGmzZtQq9evTB27FjMmTPH7es8cOAAkpKScOrUKXz99deIj493y3q4905EgDnlvgzArys/NTMUwBlVPWHCct2qc+fOWLt2LcaOHYspU6bgd7/7HcrKytyyri1btuDWW29FcXExvv76a9xyyy1uWc8VV/beZ86c6db1EJF9Gfko5GIAGwHcJCJ5IvK4iDwpIk9WDlkB4CCAXAB/B/Bbt6U1WevWrfH5559j2rRpmD17NkaNGoWTJ0+auo6vvvoKycnJaNmyJdavX4+4uDhTl1+b8PBw/OY3v8HSpUtt8ckgIrKAkTOd3PHjzjNUG2P+/PnaokUL7dixo6anpzd5eeXl5frXv/5VAwMDNTY2Vk+cOGFCSuOys7MVgM6ZM8ej6yUi94I/n6HaGI899hi2bNmCDh06YNSoUZg6dSpOnz7dqGUdOXIEd955J1544QVMmDABa9asQUREhMmJr2/w4MEYOHAgPvroI4+ul4jsgeVeTUxMDLZs2YInnngCf/vb39C7d2+8//77hufiz507h1mzZmHgwIHIzMzE/Pnz8emnnyI0NNTNyWv3yCOPYPPmzTywSuSHWO41tGzZEu+99x6ysrLQv39/PPXUU+jZsydeeOEFZGdnX/P9NGVlZdi5cyemT5+O7t27Y9q0aYiLi0NOTg4ee+wxS6+W9PDDDyMwMJB770R+SGqWlac4nU7NzMy0ZN1GqSq+/PJL/P3vf0d6ejrKysrQrl07dOrUCR07dkRJSQl27NiBixcvQkRwzz334MUXX3TbxxwbY9y4ccjJycHhw4cRGBhodRwiaiIRyVJVZ73jWO7GFBYW4rPPPsO2bduQn5+P/Px8BAUFYfDgwXA4HLjtttvQo0cPq2Ne41//+hcmTZqE9PR0jBw50uo4RNRELHcCAFy8eBGRkZEYO3YsFi5caHUcImoio+XOOXcfFxISgkmTJuGLL77AxYsXrY5DRB7CcvcD99xzD86fP4/vvvvO6ihE5CEsdz8wfPhwtGzZEsuXL7c6ChF5CMvdD4SEhGDEiBFYvny5LS81SETmY7n7iXHjxuHQoUPYvXu31VGIyANY7n5izJgxAMCpGSI/wXL3E126dEFsbCzLnchPsNz9yLhx47BhwwYUFhZaHYWI3Izl7kfGjRsHl8vl0UsLEpE1WO5+JC4uDp06deLUDJEfYLn7kYCAAIwZMwarVq1CeXm51XGIyI1Y7n5m5MiROH36NHJycqyOQkRuxHL3M8nJyQCANWvWWJqDiNyL5e5nIiMjcdNNN+Hbb7+1OgoRuRHL3Q8lJydj7dq1hi8fSETeh+Xuh5KTk3H27FnOuxP5MJa7H0pKSgLAeXciX8Zy90OcdyfyfSx3PzV8+HDOuxP5MJa7n0pOTkZxcTGys7OtjkJEbsBy91OcdyfybSx3PxUREYG+ffty3p3IR7Hc/Rjn3Yl8l6FyF5EUEdkrIrki8vtanu8uIt+KSLaIbBeRMeZHJbPdfvvtOHfuHHbs2GF1FCIyWb3lLiKBAOYAGA2gP4AHRaR/jWF/BLBUVR0AHgDwrtlByXxDhw4FAGzcuNHiJERkNiN77kMA5KrqQVUtBbAEwIQaYxRAm8rbYQCOmxeR3CUqKgqdOnXCpk2brI5CRCYzUu5dABytdj+v8rHq/gTgVyKSB2AFgKdrW5CIpIpIpohkFhQUNCIumUlEkJCQwD13Ih9kpNyllse0xv0HASxQ1a4AxgD4WESuWbaqzlNVp6o6w8PDG56WTDd06FDk5uaCb7ZEvsVIuecB6FbtfldcO+3yOIClAKCqGwGEAOhgRkByr4SEBABARkaGxUmIyExGyn0LgN4i0kNEmqHigOmyGmOOALgTAESkHyrKnbuCXsDpdCIoKIhTM0Q+pt5yV9UyAFMApAHYg4pPxewSkddEZHzlsP8E8B8isg3AYgCPqmrNqRuyoZYtW2LQoEEsdyIfE2RkkKquQMWB0uqPzah2ezeA28yNRp4ydOhQLFiwAGVlZQgKMvS/BBHZHM9QJSQkJOD8+fPYtWuX1VGIyCQsd6o6qMqpGSLfwXIn9OjRA+Hh4Sx3Ih/Ccqeqk5l4piqR72C5E4CKqZl9+/ahsLDQ6ihEZAKWOwEA4uPjAQCZmZkWJyEiM7DcCQAQGxsLgOVO5CtY7gQACAsLQ58+fVjuRD6C5U5VnE4ny53IR7DcqYrT6UReXh5OnjxpdRQiaiKWO1VxOp0AgKysLIuTEFFTsdypisPhgIhwaobIB7DcqUrr1q3Rr18/ljuRD2C501WuHFTlNzYTeTeWO13F6XTi5MmTOH6c1zgn8mYsd7rKlYOqnJoh8m4sd7rKoEGDEBgYyHIn8nIsd7pKy5YtERMTw3In8nIsd7oGD6oSeT+WO13D6XTi559/xuHDh62OQkSNxHKna8TFxQEAtm7danESImosljtdY8CAAQgMDER2drbVUYiokVjudI0WLVqgX79+3HMn8mIsd6qVw+HgnjuRF2O5U60cDgdOnDjBr/8l8lIsd6rVlcvuce+dyDux3KlWgwcPBsByJ/JWLHeqVVhYGKKjo1nuRF7KULmLSIqI7BWRXBH5fR1jJonIbhHZJSKLzI1JVnA4HPzEDJGXqrfcRSQQwBwAowH0B/CgiPSvMaY3gJcA3KaqMQCmuiEreZjD4cDBgwdx5swZq6MQUQMZ2XMfAiBXVQ+qaimAJQAm1BjzHwDmqOopAFDVfHNjkhWuHFTNycmxOAkRNZSRcu8C4Gi1+3mVj1XXB0AfEVkvIptEJMWsgGQdh8MBgF9DQOSNggyMkVoeq/l1gUEAegNIBtAVwFoRuVlVT1+1IJFUAKkA0L179waHJc/q1KkTIiMjeVCVyAsZ2XPPA9Ct2v2uAGpegy0PwJeqellVfwSwFxVlfxVVnaeqTlV1hoeHNzYzeVBsbCzLncgLGSn3LQB6i0gPEWkG4AEAy2qM+QLAcAAQkQ6omKY5aGZQsobD4cCePXtQUlJidRQiaoB6y11VywBMAZAGYA+Apaq6S0ReE5HxlcPSABSKyG4A3wJ4XlUL3RWaPCc2Nhbl5eXYsWOH1VGIqAGMzLlDVVcAWFHjsRnVbiuAaZU/5EOuHFTNzs7GkCFDLE5DREbxDFW6rhtvvBFhYWH8OCSRl2G503WJCAYPHsxyJ/IyLHeql8PhwPbt21FeXm51FCIyiOVO9Ro8eDAuXLiA/fv3Wx2FiAxiuVO9rnz9L6dmiLwHy53q1a9fPzRr1ownMxF5EZY71atZs2aIiYnhnjuRF2G5kyFXLphdcUoDEdkdy50MGTx4MAoKCnDixAmroxCRASx3MoQHVYm8C8udDBk0aBAAXjCbyFuw3MmQNm3aIDo6mnvuRF6C5U6G8WsIiLwHy50MczgcyM3NxdmzZ62OQkT1YLmTYVcOqm7fvt3iJERUH5Y7GVb9u92JyN5Y7mRYZGQkwsPDWe5EXoDlToaJSNWZqkRkbyx3ahCHw4Fdu3ahtLTU6ihEdB0sd2oQh8OBy5cvY9euXVZHIaLrYLlTg/CgKpF3YLlTg/Tq1QutW7dmuRPZHMudGiQgIACDBg1iuRPZHMudGszhcGDbtm1wuVxWRyGiOrDcqcEcDgfOnTuH3Nxcq6MQUR1Y7tRgVw6q8kvEiOyL5U4NFhMTg+DgYM67E9kYy50arFmzZujfvz/LncjGWO7UKLxgNpG9GSp3EUkRkb0ikisiv7/OuF+IiIqI07yIZEcOhwP5+fm8YDaRTdVb7iISCGAOgNEA+gN4UET61zIuFMDvAGSYHZLsh2eqEtmbkT33IQByVfWgqpYCWAJgQi3j/gzgrwAumpiPbOrKhTu2bt1qcRIiqo2Rcu8C4Gi1+3mVj1UREQeAbqq63MRsZGOhoaHo06cPsrKyrI5CRLUwUu5Sy2NVR9FEJADALAD/We+CRFJFJFNEMgsKCoynJFtyOp0sdyKbMlLueQC6VbvfFcDxavdDAdwMYI2IHAIwFMCy2g6qquo8VXWqqjM8PLzxqckW4uLikJeXh/z8fKujEFENRsp9C4DeItJDRJoBeADAsitPquoZVe2gqlGqGgVgE4DxqprplsRkG3FxcQDAvXciG6q33FW1DMAUAGkA9gBYqqq7ROQ1ERnv7oBkX1c+McNyJ7KfICODVHUFgBU1HptRx9jkpscib9CmTRseVCWyKZ6hSk0SFxfHcieyIZY7NUlcXByOHj3Kg6pENsNypybhQVUie2K5U5PwoCqRPbHcqUnCwsLQu3dvljuRzbDcqcl4UJXIflju1GRXDqryKyWI7IPlTk3Gg6pE9sNypyaLjY0FAGRm8hsniOyC5U5NFhYWhr59+2Lz5s1WRyGiSix3MkV8fDwyMjJ4TVUim2C5kyni4+ORn5+PQ4cOWR2FiMByJ5PEx8cDADIyeAldIjtguZMpBgwYgJCQEJY7kU2w3MkUwcHBcDqd2LRpk9VRiAgsdzJRfHw8srOzUVpaanUUIr/HcifTxMfH49KlS9i2bZvVUYj8HsudTMODqkT2wXIn03Tr1g2RkZGcdyeyAZY7mUZEqk5mIiJrsdzJVPHx8cjNzUVhYaHVUYj8GsudTMV5dyJ7YLmTqZxOJwICAjjvTmQxljuZKjQ0FAMHDsT69eutjkLk11juZLrExERs3LiRJzMRWYjlTqZLTExESUkJr8xEZCGWO5nu9ttvBwB8//33Fich8l8sdzJdx44d0bdvX6xdu9bqKER+i+VObpGYmIh169ahvLzc6ihEfslQuYtIiojsFZFcEfl9Lc9PE5HdIrJdRP6fiNxoflTyJomJiThz5gx27NhhdRQiv1RvuYtIIIA5AEYD6A/gQRHpX2NYNgCnqg4E8CmAv5odlLxLYmIiAM67E1nFyJ77EAC5qnpQVUsBLAEwofoAVf1WVS9U3t0EoKu5McnbdOvWDVFRUSx3IosYKfcuAI5Wu59X+VhdHgewsrYnRCRVRDJFJLOgoMB4SvJKiYmJ+P7776GqVkch8jtGyl1qeazWf60i8isATgD/XdvzqjpPVZ2q6gwPDzeekrxSYmIiCgoKsHfvXqujEPkdI+WeB6BbtftdARyvOUhERgD4A4DxqnrJnHjkzTjvTmQdI+W+BUBvEekhIs0APABgWfUBIuIAMBcVxZ5vfkzyRr169UJERATWrFljdRQiv1NvuatqGYApANIA7AGwVFV3ichrIjK+cth/A2gN4F8ikiMiy+pYHPkREcHIkSOxevVquFwuq+MQ+ZUgI4NUdQWAFTUem1Ht9giTc5GPSElJwccff4zMzEwMGTLE6jhEfoNnqJJb3XXXXRARrFq1yuooRH6F5U5u1aFDB9xyyy1YubLWT8cSkZuw3MntRo8ejc2bN/O6qkQexHInt0tJSYHL5cLq1autjkLkN1ju5Ha33HIL2rdvz3l3Ig9iuZPbBQYG4q677sKqVav4kUgiD2G5k0ekpKTgp59+wrZt26yOQuQXWO7kEaNGjQIAfmqGyENY7uQRERERiI2NxbJlPHmZyBNY7uQx999/PzIyMnDw4EGroxD5PJY7ecwDDzwAAFiyZInFSYh8H8udPCYqKgq33norFi9ebHUUIp/HciePevDBB7Fz505eOJvIzVju5FH3338/AgMDufdO5GYsd/KoTp064c4778SSJUt4bVUiN2K5k8c99NBD+PHHH5GRkWF1FCKfxXInj5s4cSKaN2+ORYsWWR2FyGcZuhITkZnatGmD8ePHY+HChfjLX/6C1q1bW5rH5XJhz5492LhxI06dOoUWLVqgRYsW6NmzJ2699VY0b97c0nxmO378eNVXMJ86dQolJSWIiopCnz59cNNNN6Ft27ZWRzSdqmLHjh04cOAA8vPzkZ+fj+bNm6NHjx6IiopC//790apVK6tjmktVLfmJi4tT8l8bN25UAPrOO+9YlmHDhg163333aVhYmAKo9adVq1Y6fvx4XbBggZaWllqWtal27dqlL7zwgg4cOLDObQWgIqJJSUk6d+5cLSwstDp2k5SXl+vKlSs1NTVVu3Tpct3tDgkJ0YkTJ+o///lPPXv2rNXRrwtAphroWJY7Web222/X7t27e7w0V69erYmJiQpA27dvr6mpqbpgwQLdt2+fFhcXa35+vh46dEiXLVumTz31lEZFRSkAvfHGG/X999/XixcvejRvU6xfv17vvvtuBaDBwcE6fPhwfeONNzQjI0MPHz6sxcXFeunSJf3hhx902bJl+sorr+hNN91UNf7JJ5/UvLw8qzejQUpLS/Xjjz/WmJgYBaChoaF633336fz583Xr1q167NgxvXTpkp4+fVq3bdumn3/+uT799NPauXNnBaBhYWH6yiuv6KlTp6zelFqx3Mn2li9frgD0448/9sj6ioqKdPLkyQpAu3btqrNmzdLi4uJ6f8/lcuny5cs1Pj5eAWhUVJSuWrXKA4kbb//+/Tp69GgFoDfccIP+6U9/0oKCAkO/63K5NCsrS5944gkNCgrSkJAQnTZtmhYVFbk5ddN9+eWX2rNnTwWgN998sy5cuFAvXbpk6HfLy8t13bp1eu+99yoAbdOmjb766qtaUlLi5tQNw3In2ysvL9eYmBgdMGCAulwut67rq6++0s6dO2tgYKD+8Y9/bNTet8vl0rS0NO3bt68C0MmTJxsuTE+5cOGCzpgxQ5s3b66hoaH65ptv6rlz5xq9vIMHD+ojjzyiAQEB2qlTJ/3kk0/c/lo1xoEDB3TcuHEKQGNiYnTZsmVaXl7e6OXl5OToxIkTFYD26tVLV69ebWLapmG5k1dYsGCBAtCVK1e6ZfllZWX60ksvVe3JZWZmNnmZJSUlOn36dA0ODtbw8HBdvny5CUmbLiMjo2pK5aGHHtLjx4+btuytW7dqXFycAtBx48bp0aNHTVt2U5SXl+vbb7+tISEh2rp1a33rrbdMneZbvXq19urVSwHoww8/bIvjECx38gqXLl3SLl266LBhw0zfIywqKtKUlBQFoKmpqabPlW/fvr3qAOXTTz9t2Z/vpaWlOn36dA0MDNSuXbtqenq6W9Zz+fJlfeutt7Rly5batm1bXbRokVvWY9SBAweqjp3cfffdbjs2cOXNPCgoSCMjI3XFihVuWY9RLHfyGvPmzVMA+t5775m2zJ07d2p0dLQGBwfr3LlzTVtuTSUlJTp16tSqvwx27NjhtnXVZvfu3VV71L/+9a89chBw//79mpCQoAD0l7/8pcf3Zl0ul86dO1dbtWqlbdq00QULFnhkqigrK6vqIG1qaqpln6phuZPXcLlcOnLkSG3VqpUeOHCgycv797//raGhoRoREaEbNmwwIWH9VqxYoR07dtSQkBCdM2eO28vG5XLp7NmzNSQkRG+44Qb97LPP3Lq+mi5fvqwzZ87UoKAgjYiI8NjU1PHjx6sOFN9555165MgRj6z3iosXL+rzzz+vIqI9evTQ7777zqPrV2W5k5c5cuSItmnTRpOSkhp9IMzlcukbb7yhIqJxcXEenxc+efJkVfHcfffdps55V3fw4EEdOXKkAtDRo0e7bT1GbN26VQcMGKAA9LHHHtPTp0+7ZT0ul0sXLFig7du31xYtWujs2bObdMC0qdatW6fR0dEqIjp16tQmHbRuKJY7eZ358+c3+sSmwsJCnTBhQtVUwfnz592QsH4ul0vfeecdDQkJ0bZt2+oHH3xg2l58WVmZvv3229qyZUtt3bq1vvvuu7b45MrFixf1pZde0oCAAI2IiNCFCxeammvfvn16xx13KABNSEjQH374wbRlN8W5c+f0t7/9rQLQbt266eeff+6R18PUcgeQAmAvgFwAv6/l+eYAPql8PgNAVH3LZLlTTS6XS8eNG6dBQUENmidfv369duvWTYODg3XWrFm2KLx9+/ZpUlKSAtDk5OQmfUrH5XLpypUr1eFwKAAdO3asx6cjjNi8ebM6nU4FoImJibp169YmLS8/P1+fe+45bd68uYaFhel7771n6d56XdatW1f118vYsWM1JyfHreszrdwBBAI4AKAngGYAtgHoX2PMbwG8X3n7AQCf1LdcljvV5vTp01WfcJk6daqWlZXVOfbo0aOampqqgYGB2rNnT92yZYsHk9avvLxc582bp+3atauaqsnKymrQ769evVqHDRtWdfLUkiVLbPHmVZeysjKdO3eutm/fXgHoqFGj9JtvvmlQ5ry8PH355Ze1VatWGhAQoJMnT7Z06smI0tJSffPNNzU0NFQB6Pjx4zUjI8Mt6zKz3BMApFW7/xKAl2qMSQOQUHk7CMDPAOR6y2W5U10uX76szzzzjALQpKQk/fDDD/XYsWOqWvHxxm+//VafffZZbd68uQYHB+vTTz/ttrleM5w5c0b//Oc/a9u2bRWADhgwQKdPn65btmy56gxZl8ulJ0+e1DVr1uizzz5bdTp8ZGSkzpkzx/CZlnZQVFSkM2fO1I4dOyoA7du3r06bNk1Xr16txcXFV5V9cXGx5uTk6Jw5czQxMVFFpGp6bffu3RZuRcMVFRXpq6++WvWG3qdPH33++ed17dq1+vPPP5vyxmxmuf8CwAfV7k8G8H9qjNkJoGu1+wcAdLjeclnuVJ+5c+dqp06dqr7cKTw8vOp2QECAPvroo/rjjz9aHdOw06dP66xZszQxMVEDAgKu+nKyqKgobdWqVdVjwcHBOn78eF20aJFeuHDB6uiNVlJSonPnztURI0Zos2bNqravefPmGhkZqREREVd9gVe/fv30tdde03379lkdvUnOnj2r7777ro4cOVKDgoKu+oKy6OhoXbx4caOXbbTcpWJs3UTkfgCjVPU3lfcnAxiiqk9XG7Orckxe5f0DlWMKaywrFUAqAHTv3j3u8OHD1103kapi+/btSE9Px+7du9GvXz8MGjQIDocDHTt2tDpeoxUUFODrr7/G0aNHcfLkSeTn5yM8PBzR0dHo2bMnEhIS0K5dO6tjmur8+fNYs2YNdu7ciaKiIhQWFkJV0atXL/Tq1Qs333wz+vbtCxGxOqqpTp8+jTVr1uDQoUPIy8vDsWPH8Pjjj2PEiBGNWp6IZKmqs95xBso9AcCfVHVU5f2XAEBV/6vamLTKMRtFJAjASQDhep2FO51OzczMNLQxRERUwWi5G7kS0xYAvUWkh4g0Q8UB02U1xiwD8Ejl7V8A+OZ6xU5ERO5V75WYVLVMRKag4qBpIID5qrpLRF5DxdzPMgD/A+BjEckFUISKNwAiIrKIocvsqeoKACtqPDaj2u2LAO43NxoRETUWL5BNROSDWO5ERD6I5U5E5INY7kREPojlTkTkg+o9icltKxYpANDYU1Q7oOL7a7yZt2+Dt+cHvH8bmN96VmzCAcGGAAADpklEQVTDjaoaXt8gy8q9KUQk08gZWnbm7dvg7fkB798G5reenbeB0zJERD6I5U5E5IO8tdznWR3ABN6+Dd6eH/D+bWB+69l2G7xyzp2IiK7PW/fciYjoOmxd7iKSIiJ7RSRXRH5fy/PNReSTyuczRCTK8ynrZiD/oyJSICI5lT+/sSJnXURkvojki8jOOp4XEflb5fZtF5FYT2esj4FtSBaRM9Vegxm1jbOKiHQTkW9FZI+I7BKRZ2oZY9vXwWB+u78GISKyWUS2VW7Dq7WMsV8XGblckxU/cNOFuW2W/1HUuGShnX4AJAKIBbCzjufHAFgJQAAMBZBhdeZGbEMygOVW57xO/kgAsZW3QwHsq+X/I9u+Dgbz2/01EACtK28HA8gAMLTGGNt1kZ333IcAyFXVg6paCmAJgAk1xkwA8FHl7U8B3Cn2uUaXkfy2pqrfo+L7+esyAcA/tMImAG1FJNIz6YwxsA22pqonVHVr5e1iAHsAdKkxzLavg8H8tlb53/Vc5d3gyp+aBytt10V2LvcuAI5Wu5+Ha/+nqBqjqmUAzgC4wSPp6mckPwDcV/mn9Kci0s0z0UxjdBvtLqHyT+6VIhJjdZi6VP6p70DFnmN1XvE6XCc/YPPXQEQCRSQHQD6A1apa52tgly6yc7nX9q5X893SyBirGMn2bwBRqjoQwNf433d+b2Hn//5GbUXF6dyDAMwG8IXFeWolIq0BfAZgqqqerfl0Lb9iq9ehnvy2fw1UtVxVBwPoCmCIiNxcY4jtXgM7l3segOp7sl0BHK9rTOWFucNgnz/B682vqoWqeqny7t8BxHkom1mMvEa2pqpnr/zJrRVXHAsWkQ4Wx7qKiASjohj/qar/t5Yhtn4d6svvDa/BFap6GsAaACk1nrJdF9m53L39wtz15q8xLzoeFfOR3mQZgF9XflpjKIAzqnrC6lANISIRV+ZGRWQIKv5NFFqb6n9VZvsfAHtU9e06htn2dTCS3wteg3ARaVt5uwWAEQB+qDHMdl1k6BqqVlAvvzC3wfy/E5HxAMpQkf9RywLXQkQWo+KTDB1EJA/AK6g4mARVfR8V19UdAyAXwAUAj1mTtG4GtuEXAJ4SkTIAJQAesPofZQ23AZgMYEflnC8AvAygO+AVr4OR/HZ/DSIBfCQigah441mqqsvt3kU8Q5WIyAfZeVqGiIgaieVOROSDWO5ERD6I5U5E5INY7kREPojlTkTkg1juREQ+iOVOROSD/j975Qt3AUQKeQAAAABJRU5ErkJggg==\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "xx = x.level()\n", "def H(w): return xx[0] + 2*sum([xx[i]*cos(i*w) for i in range(1,len(xx))])\n", "w = np.linspace(0, pi, 100)\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "plt.plot(w, [H(wi) for wi in w], 'k')\n", "plt.show()\n" ] }, { "cell_type": "markdown", "metadata": { "collapsed": true }, "source": [ "\"Creative
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). \n", "\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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 }