design-under-uncertainty.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Design Under Uncertainty\n",
    "\n",
    "This notebook was automatically generated from the Algorithms for Optimization source code. Each cell generates a figure from the original text. While this code is not optimized for use in lectures, we provide it here to be adapted for such projects. We hope you find it useful."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "include(\"support_code.jl\");"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "\tusing Distributions\n",
    "\tusing Optim\n",
    "\timport QuadGK: quadgk\n",
    "\n",
    "\tp = let\n",
    "\n",
    "\t\ta = 1\n",
    "\t\tb = 2\n",
    "\t\tc = 4.5\n",
    "\t\tf = x -> 2.5 - 2exp(-100(x-a)^2) - exp(-(x-b)^2)  - 1.4exp(-0.15(x-c)^2)\n",
    "\n",
    "\t\tp = Plots.Plot[]\n",
    "\t\tx_arr = append!(collect(range(0.0,stop=1.5,length=300)), collect(range(1.5,stop=9,length=200)[2:end]))\n",
    "\t\tpush!(p, Plots.Linear(x_arr, f.(x_arr), style=\"solid, thick, black, mark=none\", legendentry=\"noise-free\"))\n",
    "\t\tlabels=[\"very low\", \"low\", \"high\", \"very high\"]\n",
    "\t\tfor (i,σ) in enumerate([0.1,0.5,1.0,2.0])\n",
    "\t\t    P = Normal(0,σ)\n",
    "\t\t    g = x -> quadgk(z->f(x+z)*pdf(P, z), -5, 5)[1]\n",
    "\t\t    push!(p, Plots.Linear(g, (0, 9), xbins=401, style=\"solid, pastelBlue!$(100-25(i-1)), mark=none\", legendentry=labels[i]*\" noise\"))\n",
    "\t\tend\n",
    "\n",
    "\t\tmin_a = optimize(f, a-0.5, a+0.5).minimizer\n",
    "\t\tmin_b = optimize(f, b-0.5, b+0.5).minimizer\n",
    "\t\tmin_c = optimize(f, c-1.0, c+1.0).minimizer\n",
    "\n",
    "\t\tAxis(p, width=\"12cm\", height=\"6cm\", style=\"axis lines=left, xtick={$(min_a), $(min_b), $(min_c)}, xticklabels={\\$a\\$, \\$b\\$, \\$c\\$},  xticklabel style={text height=2ex}, ytick=\\\\empty, legend style={draw=none, at={(0.5,-0.25)}, anchor=north, legend columns=-1, /tikz/every even column/.append style={column sep=0.5cm}}\",\n",
    "\t\t    xlabel=L\"x\", ylabel=\"expected value\", xmin=-0.5, ymin=-0.25)\n",
    "    end\n",
    "    plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Optim\n",
    "p = let\n",
    "\tf = x -> -x*(x < 0) + x^2*(x > 0)\n",
    "\tdom = (-1,1)\n",
    "\tp = Plots.Plot[]\n",
    "\tpush!(p, Plots.Linear(f, dom, style=\"solid, black, mark=none\", legendentry=\"true\"))\n",
    "\n",
    "\tF_lin(x,ϵ) = -x + ϵ    # valid for x - ϵ ≤ 0\n",
    "\tF_qua(x,ϵ) = (x + ϵ)^2 # valid for x + ϵ > 0\n",
    "\n",
    "\tfor (i,ϵ) in enumerate([0.2,0.4,0.6,0.8,1.0,1.2,1.4,1.6,1.8,2.0])\n",
    "\t\tx_arr = Float64[dom[1]]\n",
    "\t\ty_arr = Float64[F_lin(x_arr[1], ϵ)]\n",
    "\t\tx_vertex = -(2ϵ+1)/2 + 0.5*sqrt((2ϵ+1)^2 - 4*(ϵ^2-ϵ))\n",
    "\t\tfor x in range(x_vertex, stop=dom[2], length=101)\n",
    "\t\t\tpush!(x_arr, x)\n",
    "\t\t\tpush!(y_arr, F_qua(x, ϵ))\n",
    "\t\tend\n",
    "\t\tpush!(p, Plots.Linear(x_arr, y_arr, style=\"solid, pastelBlue!$(100-9(i-1)), mark=none\", legendentry=\"\\$\\\\epsilon = $ϵ\\$\"))\n",
    "\tend\n",
    "\tAxis(p, xmin=-1, xmax=1, ymax=4, xlabel=L\"x\", ylabel=L\"y\", style=\"width=8cm, legend pos=outer north east\")\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = let\n",
    "\tZ = ϵ -> [-ϵ, ϵ]\n",
    "\tf = (x, z) -> (x+z)^2 + 6exp(-(x+z)^2)\n",
    "\n",
    "\tg = (x,ϵ) -> maximum(f(x,z) for z in range(Z(ϵ)[1], stop=Z(ϵ)[2], length=51))\n",
    "\tdom = (-3,3)\n",
    "\n",
    "\tG = GroupPlot(1,2, groupStyle=\"vertical sep=0.5cm, xlabels at=edge bottom, xticklabels at=edge bottom\", style=\"xlabel=\\$x\\$, ylabel=\\$y\\$, xmin=$(dom[1]), xmax=$(dom[2]), ymin=2, ymax=10, width=9.5cm, height=4cm\")\n",
    "\n",
    "\tp = Plots.Plot[]\n",
    "\tfor (i,ϵ) in enumerate(range(0.0,stop=3.0,length=11))\n",
    "\t\tif i == 1\n",
    "\t\t\tpush!(p, Plots.Linear(x->g(x,ϵ), dom, style=\"solid, black, mark=none\", legendentry=\"noise-free\"))\n",
    "\t\telse\n",
    "\t\t\tpush!(p, Plots.Linear(x->g(x,ϵ), dom, xbins=101, style=\"solid, pastelBlue!$(100-9(i-2)), mark=none\", legendentry=\"\\$\\\\epsilon = $(ϵ)\\$\"))\n",
    "\t\tend\n",
    "\tend\n",
    "\tpush!(p, Plots.Linear([ 2,  2],[2,10], style=\"solid, gray, mark=none, forget plot\"))\n",
    "\tpush!(p, Plots.Linear([-2, -2],[2,10], style=\"solid, gray, mark=none, forget plot\"))\n",
    "\tpush!(p, Plots.Node(L\"\\mathcal{Z}\",0,8.1))\n",
    "\tpush!(p, Plots.Node(L\"x^*\",0,5))\n",
    "\tpush!(p, Plots.Linear([-2, 2],[7.5,7.5], style=\"solid, thick, red, mark=none, forget plot\"))\n",
    "\tpush!(p, Plots.Linear([-2,-2],[7.25,7.75], style=\"solid, thick, red, mark=none, forget plot\"))\n",
    "\tpush!(p, Plots.Linear([ 2, 2],[7.25,7.75], style=\"solid, thick, red, mark=none, forget plot\"))\n",
    "\tpush!(p, Plots.Scatter([0],[f(0,0)], style=\"solid, black, mark=*, mark size=1.5, mark options={solid, draw=black, fill=black}, forget plot\"))\n",
    "\tpush!(G, Axis(p, style=\"legend pos=outer north east, legend style={draw=none, fill=none}\"))\n",
    "\n",
    "\tp = Plots.Plot[]\n",
    "\tfor (i,ϵ) in enumerate(range(0.0,stop=3.0,length=11))\n",
    "\t\tif i == 1\n",
    "\t\t\tpush!(p, Plots.Linear(x->g(x,ϵ), dom, style=\"solid, black, mark=none\"))\n",
    "\t\telse\n",
    "\t\t\tpush!(p, Plots.Linear(x->g(x,ϵ), dom, xbins=201, style=\"solid, pastelBlue!$(100-9(i-2)), mark=none\"))\n",
    "\t\tend\n",
    "\tend\n",
    "\tpush!(p, Plots.Linear([ 2,  2],[2,10], style=\"solid, gray, mark=none\"))\n",
    "\tpush!(p, Plots.Linear([-2, -2],[2,10], style=\"solid, gray, mark=none\"))\n",
    "\tpush!(p, Plots.Linear([dom[1], dom[2]],[5,5], style=\"solid, gray, mark=none\"))\n",
    "\n",
    "\tz_hi = -0.478808\n",
    "\tz_lo = -2\n",
    "\txstar = (z_hi + z_lo)/2\n",
    "\n",
    "\tpush!(p, Plots.Node(L\"\\mathcal{Z}\",-1.25,8))\n",
    "\tpush!(p, Plots.Node(L\"x^*\",xstar+0.5,f(xstar,0)+0.2))\n",
    "\tpush!(p, Plots.Linear([z_lo,z_hi],[7.5,7.5], style=\"solid, thick, red, mark=none\"))\n",
    "\tpush!(p, Plots.Linear([z_lo,z_lo],[7.25,7.75], style=\"solid, thick, red, mark=none\"))\n",
    "\tpush!(p, Plots.Linear([z_hi,z_hi],[7.25,7.75], style=\"solid, thick, red, mark=none\"))\n",
    "\tpush!(p, Plots.Scatter([xstar],[f(xstar,0)], style=\"solid, black, mark=*, mark size=1.5, mark options={solid, draw=black, fill=black}\"))\n",
    "\n",
    "\tpush!(G, Axis(p))\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Distributions\n",
    "p  = let\n",
    "\tf = x -> sin(2x)/x\n",
    "\tdom = (-10,10)\n",
    "\tp = Plots.Plot[]\n",
    "\tpush!(p, Plots.Linear(f, dom, xbins=151, style=\"solid, black, mark=none\", legendentry=\"\\$\\\\nu = 0\\$\"))\n",
    "\tfor (i,ϵ) in enumerate([0.5,1,1.5,2])\n",
    "\t    P = Normal(0, ϵ)\n",
    "\t    F = x -> quadgk(z->pdf(P, z)*f(x+z), -10, 10)[1]\n",
    "\t    push!(p, Plots.Linear(F, dom, style=\"smooth,solid, pastelBlue!$(100-10(i-1)), mark=none\", legendentry=\"\\$\\\\nu = $ϵ\\$\"))\n",
    "\tend\n",
    "\tAxis(p, xmin=dom[1], xmax=dom[2], xlabel=L\"x\", ylabel=L\"y\", style=\"width=9cm, height=4cm, legend cell align=left, legend pos=outer north east\")\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Distributions\n",
    "p  = let\n",
    "\tf = x -> isapprox(x, 0.0) ? 2.0 : sin(2x)/x\n",
    "\tdom = (-10,10)\n",
    "\tP = Normal(0, 0.1)\n",
    "\tp = Plots.Plot[]\n",
    "\tpush!(p, Plots.Linear(f, dom, xbins=151, style=\"solid, thick, black, mark=none\", legendentry=\"noise-free\"))\n",
    "\t# for (i,k) in enumerate([0.1, 0.5, 1.0, 1.5, 2.0])\n",
    "\tfor (i,k) in enumerate([1.0, 1.25, 1.5, 1.75, 2.0])\n",
    "\t    F = x -> quadgk(z->pdf(P, z)*sign(f(x+z))*abs(f(x+z))^k, -Inf, Inf)[1]\n",
    "\t    push!(p, Plots.Linear(F, dom, xbins=150, style=\"solid, pastelBlue!$(100-15(i-1)), mark=none\", legendentry=\"\\$k = $k\\$\"))\n",
    "\tend\n",
    "\tAxis(p, xmin=dom[1], xmax=dom[2], xlabel=L\"x\", ylabel=L\"y\", style=\"width=11cm, height=6cm, legend cell align=left, legend style={draw=none, at={(0.5,-0.15)}, anchor=north, legend columns=1},\")\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Distributions\n",
    "using Optim\n",
    "\n",
    "p = let\n",
    "\n",
    "\tf = x -> x^3 + 1\n",
    "\tp = Plots.Plot[]\n",
    "\tpush!(p, Plots.Linear(f, (0, 5), style=\"solid, black, mark=none\"))\n",
    "\n",
    "\tymin = -3\n",
    "\n",
    "\tδ = 0.49\n",
    "\tP = Normal(1.5,0.25)\n",
    "\tdom = quantile.(P, [0.5-δ, 0.5+δ])\n",
    "\ta, b = dom[1], dom[2]\n",
    "\tfa, fb = f(a), f(b)\n",
    "\tpush!(p, Plots.Linear(x->2pdf(P, x)+ymin, (a, b), style=\"solid, pastelBlue, mark=none\"))\n",
    "\tpush!(p, Plots.Linear([-0.5, a, a], [fa, fa, ymin], style=\"draw=none, mark=none, name path=A, forget plot\"))\n",
    "\tpush!(p, Plots.Linear([-0.5, b, b], [fb, fb, ymin], style=\"draw=none, mark=none, name path=B, forget plot\"))\n",
    "\tpush!(p, Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\"))\n",
    "\t_μ = (fa+fb)/2\n",
    "\tσ = optimize(σ->(quantile(Normal(_μ, σ), 0.5+δ) - fb)^2, 0.1, 10.0).minimizer\n",
    "\tQ = Normal(_μ, σ)\n",
    "\ty_arr = range(fa, stop=fb, length=101)\n",
    "\tx_arr = [3*(pdf(Q, y)-pdf(Q,fa))-0.5 for y in y_arr]\n",
    "\tpush!(p, Plots.Linear(x_arr, y_arr, style=\"solid, pastelBlue, mark=none\"))\n",
    "\tpush!(p, Plots.Node(\"robust\",0,fb, style=\"above right\"))\n",
    "\n",
    "\tP = Normal(4,0.25)\n",
    "\tdom = quantile.(P, [0.5-δ, 0.5+δ])\n",
    "\ta, b = dom[1], dom[2]\n",
    "\tfa, fb = f(a), f(b)\n",
    "\tpush!(p, Plots.Linear(x->2pdf(P, x)+ymin, (a, b), style=\"solid, pastelBlue, mark=none\"))\n",
    "\tpush!(p, Plots.Linear([-0.5, a, a], [fa, fa, -3], style=\"draw=none, mark=none, name path=A, forget plot\"))\n",
    "\tpush!(p, Plots.Linear([-0.5, b, b], [fb, fb, -3], style=\"draw=none, mark=none, name path=B, forget plot\"))\n",
    "\tpush!(p, Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\"))\n",
    "\t_μ = (fa+fb)/2\n",
    "\tσ = optimize(σ->(quantile(Normal(_μ, σ), 0.5+δ) - fb)^2, 0.1, 10.0).minimizer\n",
    "\tQ = Normal(_μ, σ)\n",
    "\ty_arr = range(fa, stop=fb, length=101)\n",
    "\tx_arr = [3*(pdf(Q, y)-pdf(Q,fa))-0.5 for y in y_arr]\n",
    "\tpush!(p, Plots.Linear(x_arr, y_arr, style=\"solid, pastelBlue, mark=none\"))\n",
    "\tpush!(p, Plots.Node(\"sensitive\",0,fb, style=\"above right\"))\n",
    "\n",
    "\tAxis(p, width=\"9cm\", xmin=-0.5, ymin=ymin, xlabel=L\"x\", ylabel=L\"y\", style=\"axis lines=left, axis on top, xtick=\\\\empty, ytick=\\\\empty\")\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = let\n",
    "\t_μ = x -> x^2 + 4/(1+abs(x))\n",
    "\t_σ = x -> sqrt(8/(1+abs(x)))\n",
    "\tAxis([\n",
    "\t    Plots.Linear(x -> _μ(x) + _σ(x), (-3,3), xbins=51, style=\"draw=none, mark=none, name path=A, forget plot\"),\n",
    "\t    Plots.Linear(x -> _μ(x) - _σ(x), (-3,3), xbins=51, style=\"draw=none, mark=none, name path=B, forget plot\"),\n",
    "\t    Plots.Command(\"\\\\addplot[pastelBlue!40] fill between[of=A and B];\\\\addlegendentry{\\$\\\\pm \\\\sigma\\$}\"),\n",
    "\t    Plots.Linear(x -> _μ(x), (-3,3), xbins=51, style=\"solid, black, mark=none\", legendentry=\"expected value\"),\n",
    "\t], xmin=-3, xmax=3, width=\"10cm\", height=\"5cm\", xlabel=L\"x\", ylabel=L\"y\", style=\"legend style={draw=none, at={(0.5,-0.25)}, anchor=north, legend columns=-1, /tikz/every even column/.append style={column sep=0.5cm}}\")\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "p = let\n",
    "\t_μ = x -> x^2 + 4/(1+abs(x))\n",
    "\t_σ = x -> sqrt(8/(1+abs(x)))\n",
    "\tp = Plots.Plot[]\n",
    "\tp_scatter = Plots.Plot[]\n",
    "\tfor α in range(0.0,stop=1.0,length=11)\n",
    "\t    f = x -> α*_μ(x) + (1-α)*_σ(x)\n",
    "\t    xopt = 0.0\n",
    "\t    for x in range(-3,stop=0,length=101)\n",
    "\t        if f(x) < f(xopt)\n",
    "\t            xopt = x\n",
    "\t        end\n",
    "\t    end\n",
    "\t    push!(p, Plots.Linear(f, (-3,3), xbins=151, style=\"solid, color=pastelBlue!$(round(Int, 80α+20))\"))\n",
    "\t    push!(p_scatter, Plots.Scatter([xopt, -xopt], f.([xopt, -xopt]), style=\"solid, mark=*, mark size=0.75, mark options={draw=black, fill=black}\"))\n",
    "\tend\n",
    "\tappend!(p, p_scatter)\n",
    "\tAxis(p, xmin=-3, xmax=3, width=\"10cm\", height=\"5cm\", xlabel=L\"x\", ylabel=L\"\\alpha \\mu + (1-\\alpha)\\sigma\", style=\"axis on top, clip marker paths=true, legend pos = outer north east\")\n",
    "end\n",
    "plot(p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "using Distributions\n",
    "p = let\n",
    "\txmax = 7.0\n",
    "\tymax = 0.4\n",
    "\n",
    "\tG = Truncated(Gamma(2.0,1.0), 0.0, xmax)\n",
    "\n",
    "\tμ_ = quadgk(x -> x*pdf(G, x), 0.0, xmax)[1]\n",
    "\n",
    "\tα = 0.1\n",
    "\tvalue_at_risk(α) = minimum(1-cdf(G, t) <= α ? t : Inf for t in range(0.0, stop=7.0, length=201))\n",
    "\tVaR = value_at_risk(α)\n",
    "\tcond_value_at_risk(α) = (1/α)*quadgk(t -> value_at_risk(1-t), 1-α, 1-0.001)[1]\n",
    "\tCVaR = cond_value_at_risk(α)\n",
    "\n",
    "\txq = quantile(G, 1-α)\n",
    "\n",
    "\tAxis([\n",
    "\t    Plots.Linear(z->pdf(G, z), (0, xmax), style=\"solid, pastelBlue, mark=none\"),\n",
    "\t    Plots.Linear(z->pdf(G, z), (xq, xmax), style=\"solid, pastelBlue, ultra thick, mark=none\"),\n",
    "\t    Plots.Linear([μ_,μ_], [0, ymax], style=\"solid, gray, mark=none\"),\n",
    "\t    Plots.Linear([VaR,VaR], [0, ymax], style=\"solid, gray, mark=none\"),\n",
    "\t    Plots.Linear([CVaR,CVaR], [0, ymax], style=\"solid, gray, mark=none\"),\n",
    "\t    Plots.Node(\"top \\$1-\\\\alpha\\$ quantile\", (xq+xmax)/2+0.35, 0.08),\n",
    "\t], xlabel=L\"y\", ylabel=L\"p(y)\", style=\"xtick={$μ_,$VaR, $CVaR}, xticklabels={expected value, VaR, CVaR}\", ymin=0, ymax=ymax, xmin=0, xmax=xmax, width=\"12cm\", height=\"4cm\"\n",
    "\t)\n",
    "end\n",
    "plot(p)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.5.1",
   "language": "julia",
   "name": "julia-1.5"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}