import%20marimo%0A%0A__generated_with%20%3D%20%220.15.2%22%0Aapp%20%3D%20marimo.App()%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22!%5BMOSEK%20ApS%5D(https%3A%2F%2Fwww.mosek.com%2Fstatic%2Fimages%2Fbranding%2Fwebgraphmoseklogocolor.png%20)%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%20Hard%20Uncertain%20Convex%20Inequalities%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22Robust%20optimization%20is%20a%20popular%20technique%20for%20dealing%20with%20uncertainty%20in%20optimization%20problems.%20Methods%20of%20reformulating%20uncertain%20problems%20into%20their%20robust%20counterparts%20exist%20for%20a%20wide%20range%20of%20problems.%20However%2C%20problems%20with%20constraints%20that%20are%20convex%20in%20both%20their%20optimization%20variables%20and%20the%20uncertain%20parameters%20are%2C%20in%20general%2C%20harder%20to%20solve.%20Roos%20et.%20al.%20(2018)%20(http%3A%2F%2Fwww.optimization-online.org%2FDB_HTML%2F2018%2F06%2F6679.html)%20provide%20a%20way%20to%20approximate%20this%20popular%20class%20of%20problems.%20Implementing%20their%20approximations%20in%20Fusion%20is%20going%20to%20be%20the%20goal%20of%20this%20notebook.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20Adjustable%20Robust%20Optimization%20problem%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20Consider%20the%20uncertain%20optimization%20problem%3A%20%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%20%5Cmathrm%7Bminimize%7D%20%5Cquad%20c%5ET%20x%20%0A%20%20%20%20%24%24%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%20%5Ctext%7Bsubject%20to%7D%20%5Cquad%20f(A(x)%20%5Czeta%20%2B%20b(x))%20%5Cleq%200%20%5Cquad%20%5Cforall%20%5Czeta%20%5Cin%20U%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24f%3A%20%5Cmathbb%7BR%7D%5Ep%20%5Cto%20%5Cmathbb%7BR%7D%24%20is%20convex%20and%20closed%2C%20%0A%20%20%20%20%24A%3A%20%5Cmathbb%7BR%7D%5En%20%5Cto%20%5Cmathbb%7BR%7D%5E%7Bp%20%5Ctimes%20L%7D%24%2C%20%0A%20%20%20%20%24b%3A%20%5Cmathbb%7BR%7D%5En%20%5Cto%20%5Cmathbb%7BR%7D%5Ep%24%20are%20affine%2C%20%0A%20%20%20%20and%20%24U%24%20is%20the%20uncertainty%20set%2C%20which%20is%20a%20non-empty%20polyhedron%20given%20by%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20U%20%3D%20%5C%7B%5Czeta%20%5Cin%20%5CBbb%20R%5EL%20%7C%20D%5Czeta%20%3D%20d%5C%7D%0A%20%20%20%20%24%24%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22for%20some%20%24D%20%5Cin%20%5CBbb%20R%5E%7Bq%5Ctimes%20L%7D%24%20and%20%24d%20%5Cin%20%5CBbb%20R%5Eq%24.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20The%20above-stated%20constraint%20will%20be%20satisfied%20if%20and%20only%20if%20%24x%24%20satisfies%20the%20following%20adjustable%20robust%20optimization%20constraints%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cforall%20w%20%5Cin%20%5Ctext%7Bdom%7D%20f%5E*%2C%20%5Cexists%20%5Clambda%20%5Cin%20%5CBbb%20R%5Eq%20%3A%20%5Cbegin%7Bcases%7D%0A%20%20%20%20d%5ET%5Clambda%20%2B%20w%5ETb(x)%20-%20f%5E*(x)%5Cleq%200%20%5C%5C%0A%20%20%20%20D%5ET%5Clambda%20%3D%20A(x)%5ETw%2C%0A%20%20%20%20%5Cend%7Bcases%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24f%5E*(w)%3A%5CBbb%20R%5Ep%20%5Cmapsto%20%5CBbb%20R%24%20is%20the%20convex%20conjugate%20of%20the%20original%20function%20%24f%24%20in%20the%20constraint.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%23%20Note%3A%20Since%20we%20have%20%24%5Czeta%20%5Cin%20%5CBbb%20R%5EL%24%2C%20we%20get%20an%20equality%20as%20the%20second%20constraint.%20If%20%24%20%5Czeta%20%24%20has%20a%20non-negativity%20constraint%20in%20%24U%24%2C%20i.e.%20%24%5Czeta%20%5Cin%20%5CBbb%20R_%2B%5EL%24%2C%20then%20the%20equality%20is%20replaced%20by%20%24D%5ET%20%5Clambda%20%5Cgeq%20A(x)%5ETw%24.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20Safe%20Approximations%3A%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22Roos%20et.%20al.%20provide%20two%20safe%20approximations%20to%20the%20ARO%20given%20above%3A%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%20Approximation%201%3A%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20If%20there%20exists%20%24u%5Cin%20%5CBbb%20R%5Eq%24%20and%20%24V%5Cin%20%5CBbb%20R%5E%7Bq%5Ctimes%20p%7D%24%20for%20a%20given%20%24x%5Cin%20%5CBbb%20R%5En%24%20such%20that%20%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Bcases%7D%0A%20%20%20%20d%5ETu%20%2B%20f%5Cbig(b(x)%2BV%5ETd%5Cbig)%20%5Cleq%200%20%5C%5C%0A%20%20%20%20%5Cdelta%5E*%5Cbig(A_i(x)-V%5ETD_i%7Cdom%5C%2Cf%5E*%5Cbig)%20%3D%20D%5ET_iu%20%5Cquad%20i%3D1%2C...%2CL.%0A%20%20%20%20%5Cend%7Bcases%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20holds%2C%20then%20%24x%24%20also%20satisfies%20the%20original%20statement%20of%20the%20problem.%20Here%2C%20%24A_i%24%20stands%20for%20the%20%24i%24-th%20row%20of%20the%20matrix%20%24A%24.%20%24%5Cdelta%5E*%24%20is%20the%20support%20function.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%20Approximation%202%3A%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20If%20there%20exists%20%24u%5Cin%20%5CBbb%20R%5Eq%24%2C%20%24V%5Cin%20%5CBbb%20R%5E%7Bq%5Ctimes%20p%7D%24%20and%20%24r%5Cin%20%5CBbb%20R%5Eq%24%20for%20a%20given%20%24x%5Cin%20%5CBbb%20R%5En%24%20such%20that%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Bcases%7D%0A%20%20%20%20d%5ETu%20%2B%20(1%2Bd%5ETr)f%5CBig(%20%5Cfrac%7BV%5ETd%2Bb(x)%7D%7B1%2Bd%5ETr%7D%5CBig)%20%5Cleq%200%20%5C%5C%0A%20%20%20%201%2Bd%5ETr%20%5Cgeq%200%20%5C%5C%0A%20%20%20%20-D_i%5ETu%20%2B%20(-D_i%5ETr)f%5CBig(%20%5Cfrac%7BA_i(x)-V%5ETD_i%7D%7B-D_i%5ETr%7D%20%5CBig)%20%3D%200%20%5Cquad%20i%3D%201%2C...%2CL%5C%5C%0A%20%20%20%20-D_i%5ETr%20%5Cgeq%200%20%0A%20%20%20%20%5Cend%7Bcases%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20holds%2C%20then%20%24x%24%20also%20satisfies%20the%20original%20problem%20statement.%20Moreover%2C%20this%20is%20a%20tighter%20approximation%20than%20the%20first%20approximation.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20Numerical%20example%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20Roos%20et.%20al.%20test%20the%20above-mentioned%20approach%20by%20solving%20the%20following%20problem%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Ctext%7BMinimize%7D%5Cquad%20c%5ETx%0A%20%20%20%20%24%24%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Ctext%7BSubject%20to%7D%5Cquad%20h%5CBigg(%0A%20%20%20%20%5Cbegin%7Bpmatrix%7D%0A%20%20%20%20%5Cbig(%20-%5Cmathbf%7B1%7D%20%2B%20B_i%5E%7B(1)%7D%5Czeta%20%5Cbig)%5ET%20x%20%5C%5C%0A%20%20%20%20%5Cbig(%20-%5Cmathbf%7B1%7D%20%2B%20B_i%5E%7B(2)%7D%5Czeta%20%5Cbig)%5ET%20x%0A%20%20%20%20%5Cend%7Bpmatrix%7D%20%5CBigg)%20%5Cleq%200%20%20%5Cquad%20%5Cquad%20%5Cforall%20%5Czeta%20%5Cin%20U%2C%5C%2C%20i%3D1%2C...%2C%20m%0A%20%20%20%20%24%24%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20where%20%24c%3D%5Cmathbf%7B1%7D%5Cin%20%5CBbb%20R%5En%24%20is%20the%20all%20ones%20vector%2C%20and%20%24B_i%5E%7B(1)%7D%2CB_i%5E%7B(2)%7D%20%5Cin%20%5CBbb%20R%5E%7Bn%5Ctimes%20L%7D%24%20are%20randomly%20generated%20sparse%20matrices%20(sparsity%20density%20%3D%200.1)%2C%20with%20non-zero%20elements%20that%20lie%20in%20the%20interval%20%24%5B0%2C1%5D%24.%20The%20uncertainty%20set%20is%20assumed%20to%20be%20a%20box%2C%20that%20is%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20U%20%3D%20%5Cbig%5C%7B%20%5Czeta%20%5Cin%20%5CBbb%20R%5EL%20%5C%2C%5Cbig%7C%5C%2C%5C%2C%20%5Cleft%5ClVert%20%5Czeta%20%5Cright%5CrVert_%5Cinfty%20%5Cleq%201%20%5Cbig%5C%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20which%20can%20be%20re-stated%20as%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20U%20%3D%20%5Cbig%5C%7B%20%5Czeta%20%5Cin%20%5Cmathbb%7BR%7D%5EL%20%5C%2C%5Cbig%7C%5C%2C%20D%5Czeta%20%5Cleq%20d%20%5Cbig%5C%7D.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20with%0A%0A%20%20%20%20%24%24%0A%20%20%20%20D%20%3D%20%5Cbegin%7Bbmatrix%7D%0A%20%20%20%20%5Cmathbf%7BI%7D%20%5C%5C%0A%20%20%20%20%5Cmathbf%7B-I%7D%0A%20%20%20%20%5Cend%7Bbmatrix%7D%20%20%0A%20%20%20%20%5Cquad%0A%20%20%20%20d%20%3D%20%5Cmathbf%7B1%7D%20%5Cin%20%5CBbb%20R%5EL%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24I%20%5Cin%20%5CBbb%20R%5E%7BL%5Ctimes%20L%7D%24%20is%20the%20identity%20matrix.%20Therefore%2C%20we%20see%20that%20%24q%3D2L%24.%20Moreover%2C%20by%20making%20the%20following%20substitutions%3A%0A%0A%20%20%20%20%24%24A(x)%20%3D%0A%20%20%20%20%5Cbegin%7Bbmatrix%7D%0A%20%20%20%20%5Cbig(%20B_i%5E%7B(1)%20T%7D%20x%20%5Cbig)%20%5ET%20%5C%5C%0A%20%20%20%20%5Cbig(%20B_i%5E%7B(2)%20T%7D%20x%20%5Cbig)%20%5ET%0A%20%20%20%20%5Cend%7Bbmatrix%7D%20%0A%20%20%20%20%5Cquad%0A%20%20%20%20b(x)%20%3D%20%5Cbegin%7Bbmatrix%7D%0A%20%20%20%20-%5Cmathbf%7B1%7D%5ET%20x%20%5C%5C%20%0A%20%20%20%20-%5Cmathbf%7B1%7D%5ET%20x%0A%20%20%20%20%5Cend%7Bbmatrix%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20we%20can%20express%20the%20example%20constraints%20in%20the%20form%3A%20%24f%5Cbig(A(x)%5Czeta%20%2B%20b(x)%5Cbig)%24%2C%20where%20the%20function%20%24f%24%20is%20now%20the%20log-sum-exp%20(%24h%24)%20function%20with%20%24p%3D2%24.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%201.)%20Exact%20solution%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20We%20can%20solve%20the%20robust%20optmization%20problem%20exactly%2C%20by%20evaluating%20the%20constraints%20at%20every%20vertex%20of%20the%20uncertainty%20box%20set.%20Therefore%2C%20we%20end%20up%20with%20%24m%5Ctimes2%5EL%24%20log-sum-exp%20inequalities.%20Each%20inequality%20can%20be%20expressed%20as%3A%20%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5CBig(h%5E%7B(j)%7D_1%2Ch%5E%7B(j)%7D_2%2Ch%5E%7B(j)%7D_3%20%5CBig)%20%5Cin%20K_%7Bexp%7D%20%5Cquad%2C%20%5Cquad%20%20%5Csum_%7Bj%7D%20h_1%5E%7B(j)%7D%20%5Cleq%201%20%5Cquad%20%5Cquad%20j%3D1%2C...%2Cp%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24h_2%5E%7Bi%7D%20%3D%201%24%20and%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20h%5E%7B(j)%7D_3%20%3D%20-%5Cmathbf%7B1%7D%5ET%20x%20%20%2B%20%5Cbig(B_i%5E%7B(j)T%7Dx%5Cbig)%5ET%20%5Czeta%20%20%5Cquad%20i%3D1%2C...%2Cm%0A%20%20%20%20%24%24%0A%0A%20%20%20%20for%20a%20given%20uncertainty%20vector%20%24%5Czeta%20%5Cin%20U%24.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20from%20numpy%20import%20random%0A%20%20%20%20from%20mosek.fusion%20import%20Model%2C%20Domain%2C%20Expr%2C%20ObjectiveSense%2C%20Matrix%2C%20Var%0A%20%20%20%20import%20mosek.fusion.pythonic%0A%20%20%20%20import%20sys%0A%20%20%20%20import%20scipy%0A%20%20%20%20from%20scipy%20import%20sparse%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20plt.rcParams%5B'figure.figsize'%5D%20%3D%20%5B16%2C%209%5D%0A%20%20%20%20return%20Domain%2C%20Expr%2C%20Matrix%2C%20Model%2C%20ObjectiveSense%2C%20np%2C%20plt%2C%20scipy%0A%0A%0A%40app.cell%0Adef%20_(scipy)%3A%0A%20%20%20%20%23Function%20that%20generates%20a%20random%20sparse%20matrix%20B%2C%20with%20the%20specified%20sparsity%20density%2C%20and%20dimensions%20n*L.%0A%20%20%20%20%23Please%20note%20that%20we%20generate%20transposed%20matrices.%0A%20%20%20%20def%20Sparse_make(n%2CL%2Cs%2Cran_stat)%3A%0A%20%20%20%20%20%20%20%20return%20scipy.sparse.random(L%2Cn%2Cdensity%20%3D%20s%2Crandom_state%20%3D%20ran_stat).todense()%0A%20%20%20%20return%20(Sparse_make%2C)%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20%23Function%20to%20set%20the%20dimensional%20constants%2C%20matrix%20D%2C%20and%20the%20vector%20d.%20%0A%20%20%20%20def%20constants(L)%3A%0A%20%20%20%20%20%20%20%20q%20%3D%202*L%0A%20%20%20%20%20%20%20%20d%20%3D%20np.ones(q)%0A%20%20%20%20%20%20%20%20D%20%3D%20np.concatenate((np.eye(L)%2C-1*np.eye(L))%2Caxis%3D0)%0A%20%20%20%20%20%20%20%20Dt%20%3D%20np.transpose(D)%0A%20%20%20%20%20%20%20%20return(q%2Cd%2CD%2CDt)%0A%20%20%20%20return%20(constants%2C)%0A%0A%0A%40app.cell%0Adef%20_(Matrix%2C%20np)%3A%0A%20%20%20%20%23Making%20the%20vertex%20vectors%20of%20the%20uncertainty%20box.%20This%20will%20be%20used%20in%20the%20exact%20solution.%0A%20%20%20%20def%20Box_vertices(L)%3A%0A%20%20%20%20%20%20%20%20U%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20for%20i%20in%20range(2**L)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20vertex%20%3D%20%5B2*int(b)%20-%201%20for%20b%20in%20list(format(i%2C'b').zfill(L))%5D%0A%20%20%20%20%20%20%20%20%20%20%20%20U.append(vertex)%0A%20%20%20%20%20%20%20%20U_matrix%20%3D%20Matrix.dense(np.transpose(np.asarray(U)))%0A%20%20%20%20%20%20%20%20return(U_matrix)%0A%20%20%20%20return%20(Box_vertices%2C)%0A%0A%0A%40app.cell%0Adef%20_(Box_vertices%2C%20Domain%2C%20Expr%2C%20Model%2C%20ObjectiveSense%2C%20constants)%3A%0A%20%20%20%20def%20Exact_sol(Bt_p%2Cn%2Cm%2CL%2Cp)%3A%20%20%20%20%0A%20%20%20%20%20%20%20%20U%20%3D%20Box_vertices(L)%0A%20%20%20%20%20%20%20%20q%2Cd%2CD%2Cdt%20%3D%20constants(L)%0A%0A%20%20%20%20%20%20%20%20with%20Model('EXACT')%20as%20M%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23M.setLogHandler(sys.stdout)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.setSolverParam('numThreads'%2C%201)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20M.variable('x'%2Cn)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Objective%0A%20%20%20%20%20%20%20%20%20%20%20%20M.objective('OBJ_exact'%2CObjectiveSense.Minimize%2CExpr.sum(x))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23b(x)%2Cdim%20%3D%20%5Bp%2C2%5EL%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20bx%20%3D%20Expr.repeat(Expr.repeat(Expr.repeat(-Expr.sum(x)%2Cp%2C0)%2C2**L%2C1)%2Cm%2C2)%20%20%20%20%20%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%5BA(x).(uncertainty%20vector)%5D%20%2C%20dim%20%3D%20%5Bp%2C2%5EL%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20Ax_U%20%3D%20Expr.stack(2%2C%5BExpr.hstack(%5Bb%20%40%20x%20for%20b%20in%20B_set%5D).T%20%40%20U%20for%20B_set%20in%20Bt_p%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20h1%20%3D%20M.variable('h1'%2C%5Bp%2C2**L%2Cm%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20h2%20%3D%20Expr.constTerm(%5Bp%2C2**L%2Cm%5D%2C1.0)%0A%20%20%20%20%20%20%20%20%20%20%20%20h3%20%3D%20Ax_U%2Bbx%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20Cone%20%3D%20M.constraint('Exp1'%2CExpr.stack(3%2C%5Bh1%2Ch2%2Ch3%5D)%2CDomain.inPExpCone())%0A%20%20%20%20%20%20%20%20%20%20%20%20Line%20%3D%20M.constraint('Lin1'%2CExpr.sum(h1%2C0)%3C%3D1.0)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20%20%20%20%20return(M.primalObjValue()%2CM.getSolverDoubleInfo('optimizerTime'))%0A%20%20%20%20return%20(Exact_sol%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%202.)%20Approximation%201%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20It%20is%20easy%20to%20model%20the%20first%20approximation%20by%20using%20the%20previously%20mentioned%20expressions%20for%20%24A(x)%24%2C%20%24b(x)%24%2C%20%24D%24%20and%20%24d%24.%20We%20see%20that%20for%20the%20first%20constraint%20in%20the%20approximation%2C%20i.e%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20d%5ETu%20%2B%20h(V%5ETd%20%2B%20b(x))%20%5Cleq%200%0A%20%20%20%20%24%24%0A%0A%20%20%20%20we%20can%20write%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Csum_%7Bk%3D1%7D%5E%7Bp%7De%5E%7B%5Cbig(%20V_k%5ETd%20%2B%20b(x)_k%20%2B%20d%5ETu%20%5Cbig)%7D%20%5Cleq%201%0A%20%20%20%20%24%24%0A%0A%20%20%20%20since%20%24h%24%20is%20a%20log-sum-exp%20function.%20Note%20that%20k%20is%20the%20column%20number%20in%20%24V%24%20(two%20dimensional%20variable)%20and%20the%20row%20number%20in%20%24b(x)%24.%20The%20above%20expression%20can%20be%20easily%20expressed%20in%20Fusion%2C%20using%20exponential%20cones%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbig(h_1%5E%7B(k)%7D%2Ch_2%5E%7B(k)%7D%2Ch_1%5E%7B(k)%7D%5Cbig)%20%5Cin%20K_%7Bexp%7D%20%5Cquad%20%5Cquad%20%5Csum_%7Bk%3D1%7D%5E%7Bp%7D%20h_1%5E%7B(k)%7D%20%5Cleq%201%20%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24h_2%5E%7B(k)%7D%20%3D%201%24%20and%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20h_3%5E%7B(k)%7D%20%3D%20V_k%5ETd%20%2B%20b(x)_k%20%2B%20d%5ETu%20%5Cquad%20%5Cquad%20k%3D1%2C...%2Cp%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20second%20constraint%20in%20the%20approximation%2C%20i.e.%20the%20equality%2C%20involves%20the%20support%20function%20of%20the%20domain%20of%20%24f%5E*%24%20(i.e.%20the%20conjugate).%20For%20the%20case%20of%20log-sum-exp%20%24(h)%24%2C%20we%20find%20that%20the%20domain%20is%20a%20simplex%20and%20the%20support%20function%20for%20a%20simplex%20is%20the%20maximum%20function.%20Hence%2C%20the%20second%20constraint%20is%20easily%20expressed%20as%20a%20set%20of%20linear%20equations%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cmax_%7Bk%7D%20%5Cbig%5C%7B%20A_%7Bki%7D(x)-V%5ET_k%20D_i%20%5Cbig%5C%7D%20%3D%20D_i%5ET%20u%20%20%5Cquad%20%5Cquad%20i%20%3D%201%2C...%2CL%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24k%20%3D%201%2C...%2Cp%24%20and%20%24A_%7Bki%7D%24%20denotes%20the%20element%20in%20%24A(x)%24%20at%20the%20%24k%24-th%20row%20and%20%24i%24-th%20column.%20%24V_k%24%20denotes%20the%20%24k%24-th%20column%20in%20%24V%24.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Domain%2C%20Expr%2C%20Matrix%2C%20Model%2C%20ObjectiveSense%2C%20constants%2C%20np)%3A%0A%20%20%20%20def%20approx_1(Bt_p%2Cn%2Cm%2CL%2Cp)%3A%0A%20%20%20%20%20%20%20%20q%2Cd%2CD%2CDt%20%3Dconstants(L)%0A%20%20%20%20%20%20%20%20d%20%3D%20np.reshape(np.array(d)%2C%20(len(d)%2C%201))%0A%0A%20%20%20%20%20%20%20%20with%20Model('APPROX_1')%20as%20M%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23M.setLogHandler(sys.stdout)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.setSolverParam('numThreads'%2C%201)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20M.variable('x'%2Cn)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Objective%0A%20%20%20%20%20%20%20%20%20%20%20%20M.objective('OBJ_approx1'%2CObjectiveSense.Minimize%2CExpr.sum(x))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20u%20%3D%20M.variable('u'%2C%5Bq%2Cm%5D%2CDomain.greaterThan(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20Vt%20%3D%20M.variable('Vt'%2C%5Bp%2Cq%2Cm%5D%2CDomain.greaterThan(0.0))%20%20%20%20%20%20%20%20%20%20%20%20%20%23Note%20that%20we%20define%20the%20transpose%20of%20V.%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23b(x)%2C%20dim%20%3D%20%5Bp%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20bx%20%3D%20Expr.repeat(Expr.repeat(-Expr.sum(x)%2Cp%2C0)%2Cm%2C1)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23A(x)%2C%20dim%20%3D%20%5Bp%2CL%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20Ax%20%3D%20Expr.stack(2%2C%5BExpr.hstack(%5BbT%5Bi%5D%40x%20for%20i%20in%20range(p)%5D).T%20for%20bT%20in%20Bt_p%5D)%20%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20h1%20%3D%20M.variable('h1'%2C%5Bp%2Cm%5D%2CDomain.greaterThan(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20h2%20%3D%20Expr.constTerm(%5Bp%2Cm%5D%2C1.0)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Matrix%20mul(V%5ET%2Cd)%2C%20dim%20%3D%20%5Bp%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20Vt_d%20%3D%20Expr.hstack(%5BVt%5B%3Ap%2C%3Aq%2Ci%5D.reshape(%5Bp%2Cq%5D)%20%40%20d%20for%20i%20in%20range(m)%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23(d%5ET%20.%20u)%2C%20dim%20%3D%20%5Bp%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20dt_u%20%3D%20Expr.repeat(d.T%40u%2Cp%2C0)%0A%20%20%20%20%20%20%20%20%20%20%20%20h3%20%3D%20Vt_d%20%2B%20bx%20%2B%20dt_u%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Exponential%20Cone%0A%20%20%20%20%20%20%20%20%20%20%20%20Cone%20%3D%20M.constraint('ExpCone'%2CExpr.stack(2%2C%5Bh1%2Ch2%2Ch3%5D)%2CDomain.inPExpCone())%20%20%20%20%20%20%20%20%20%20%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Linear%20constraint%20on%20the%20vector%20h1.%0A%20%20%20%20%20%20%20%20%20%20%20%20LinIneq%20%3D%20M.constraint('LinIneq'%2CExpr.sum(h1%2C0)%3C%3D1.0)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Matrix%20mult(V%5ET%2CD)%2C%20dim%20%3D%20%5Bp%2CL%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20Vt_D%20%3D%20Expr.stack(2%2C%5BVt%5B%3Ap%2C%3Aq%2Ci%5D.reshape(%5Bp%2Cq%5D)%20%40%20Matrix.sparse(D)%20for%20i%20in%20range(m)%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Matrix%20mult(D%5ET%2Cu)%2C%20dim%20%3D%20%5Bp%2CL%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20Dt_u%20%3D%20Expr.repeat(Expr.reshape(Dt%20%40%20u%2C%5B1%2CL%2Cm%5D)%2Cp%2C0)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Second%20constraint%2C%20i.e.%20the%20equality.%20%0A%20%20%20%20%20%20%20%20%20%20%20%20Eq%20%3D%20M.constraint('Eq'%2CAx%20%3D%3D%20Vt_D%20%2B%20Dt_u)%20%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20%20%20%20%20return(M.primalObjValue()%2CM.getSolverDoubleInfo('optimizerTime'))%0A%20%20%20%20return%20(approx_1%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%23%203.)%20Approximation%202%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20For%20the%20second%20approximation%2C%20we%20can%20manipulate%20the%20first%20constraint%20into%20the%20following%20form%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Csum_%7Bk%3D1%7D%5E%7Bp%7De%5E%7B%5CBig(%20%5Cfrac%7BV_k%5ETd%20%2B%20b(x)_k%20%2B%20d%5ETu%7D%7B1%2Bd%5ETr%7D%20%5CBig)%7D%20%5Cleq%201%0A%20%20%20%20%24%24%0A%0A%20%20%20%20which%20then%20gives%20us%20the%20following%20exponential%20cones%2C%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbig(h_%7B11%7D%5E%7B(k)%7D%2Ch_%7B12%7D%5E%7B(k)%7D%2C%20h_%7B13%7D%5E%7B(k)%7D%5Cbig)%20%5Cin%20K_%7Bexp%7D%20%5Cquad%20%5Cquad%20%5Csum_%7Bk%3D1%7D%5E%7Bp%7D%20h_%7B11%7D%5E%7B(k)%7D%20-%20(1%2Bd%5ETr)%20%5Cleq%200%20%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24h_%7B12%7D%5E%7B(k)%7D%20%3D%20%5Cbig(1%2Bd%5ETr%5Cbig)%24%20and%20%24h_%7B13%7D%5E%7B(k)%7D%20%3D%20%5Cbig(V_k%5ETd%20%2B%20b(x)_k%20%2B%20d%5ETu%5Cbig)%24.%0A%20%20%20%20Following%20the%20same%20approach%20for%20the%20third%20inequality%20gives%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbigg(h_%7B21%7D%5E%7B(k)%7D%2Ch_%7B22%7D%5E%7B(k)%7D%2Ch_%7B23%7D%5E%7B(k)%7D%20%5Cbigg)%20%5Cin%20K_%7Bexp%7D%20%5Cquad%20%5Cquad%20%5Csum_%7Bk%3D1%7D%5E%7Bp%7D%20h_%7B21%7D%5E%7B(k)%7D%20%2B%20D_i%5ETr%20%5Cleq%200%20%20%5Cquad%20%5Cquad%20i%20%3D%201%2C...%2CL%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24h_%7B22%7D%5E%7B(k)%7D%20%3D%20-D_i%5ETr%24%20and%20%24h_%7B23%7D%5E%7B(k)%7D%20%3D%20%5Cbig(%20A_i(x)%20-%20V%5ETD_i%20-D_i%5ETu%20%5Cbig)%24.%0A%0A%20%20%20%20Note%20that%20%24h_%7B12%7D%5E%7B(k)%7D%24%20and%20%24h_%7B22%7D%5E%7B(k)%7D%24%20are%20included%20in%20the%20cone%20domain%20and%20therefore%20the%20second%20and%20the%20fourth%20inequalities%20in%20the%20statement%20for%20second%20approximation%20are%20taken%20care%20of.%0A%0A%20%20%20%20%23%23%23%23%23%20Mosek%20highlight%3A%20Roos%20et.%20al.%20mention%20that%20interior%20point%20methods%20can%20have%20trouble%20solving%20the%20second%20approximation%20(very%20small%20denominator%20case).%20However%2C%20Mosek%20successfully%20solves%20every%20single%20instance%20of%20the%20problem%20without%20any%20issues.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(Domain%2C%20Expr%2C%20Matrix%2C%20Model%2C%20ObjectiveSense%2C%20constants%2C%20np)%3A%0A%20%20%20%20def%20approx_2(Bt_p%2Cn%2Cm%2CL%2Cp)%3A%0A%20%20%20%20%20%20%20%20q%2Cd%2CD%2CDt%20%3Dconstants(L)%0A%20%20%20%20%20%20%20%20d%20%3D%20np.reshape(np.array(d)%2C%20(len(d)%2C%201))%0A%0A%20%20%20%20%20%20%20%20with%20Model('APPROX_2')%20as%20M%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23M.setLogHandler(sys.stdout)%0A%20%20%20%20%20%20%20%20%20%20%20%20M.setSolverParam('numThreads'%2C%201)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20x%20%3D%20M.variable('x'%2Cn)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Objective%0A%20%20%20%20%20%20%20%20%20%20%20%20M.objective('OBJ_approx2'%2CObjectiveSense.Minimize%2CExpr.sum(x))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23b(x)%2C%20dim%20%3D%20%5Bp%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20bx%20%3D%20Expr.repeat(Expr.repeat(-Expr.sum(x)%2Cp%2C0)%2Cm%2C1)%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%23A(x)%2C%20dim%20%3D%20%5Bp%2CL%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20Ax%20%3D%20Expr.stack(2%2C%5BExpr.hstack(%5BbT%5Bi%5D%40x%20for%20i%20in%20range(p)%5D).T%20for%20bT%20in%20Bt_p%5D)%20%20%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20u%20%3D%20M.variable('u'%2C%5Bq%2Cm%5D%2CDomain.greaterThan(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20r%20%3D%20M.variable('r'%2C%5Bq%2Cm%5D%2CDomain.greaterThan(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20Vt%20%3D%20M.variable('Vt'%2C%5Bp%2Cq%2Cm%5D%2CDomain.greaterThan(0.0))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20h_11%20%3D%20M.variable('h_11'%2C%5Bp%2Cm%5D%2CDomain.greaterThan(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20%231%20%2B%20(d%5ET.r)%2C%20dim%20%3D%20%5Bm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20dt_r_1%20%3D%20Expr.constTerm(%5Bm%5D%2C1.0)%20%2B%20(d.T%20%40%20r).F%0A%20%20%20%20%20%20%20%20%20%20%20%20h_12%20%3D%20Expr.repeat(dt_r_1%2Cp%2C1).T%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Matrix%20mult(V%5ET%2Cd)%2C%20dim%20%3D%20%5Bp%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20Vt_d%20%3D%20Expr.hstack(%5BVt%5B%3Ap%2C%3Aq%2Ci%5D.reshape(%5Bp%2Cq%5D)%20%40%20d%20for%20i%20in%20range(m)%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23(d%5ET%2Cu)%2C%20dim%20%3D%20%5Bp%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20dt_u%20%3D%20Expr.repeat(d.T%20%40%20u%2Cp%2C0)%0A%20%20%20%20%20%20%20%20%20%20%20%20h_13%20%3D%20Vt_d%20%2B%20dt_u%20%2B%20bx%0A%20%20%20%20%20%20%20%20%20%20%20%20%23First%20set%20of%20cones.%0A%20%20%20%20%20%20%20%20%20%20%20%20Cone1%20%3D%20M.constraint('Cone1'%2CExpr.stack(2%2C%5Bh_11%2Ch_12%2Ch_13%5D)%2CDomain.inPExpCone())%0A%20%20%20%20%20%20%20%20%20%20%20%20LinIneq1%20%3D%20M.constraint('LinIneq1'%2CExpr.sum(h_11%2C0)%20%3C%3D%20dt_r_1)%0A%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20h_21%20%3D%20M.variable('h_21'%2C%5Bp%2CL%2Cm%5D%2CDomain.greaterThan(0.0))%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Matrix%20mult(D_T%2Cr)%2C%20dim%20%3D%20%5BL%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20Dt_r%20%3D%20Dt%20%40%20r%0A%20%20%20%20%20%20%20%20%20%20%20%20h_22%20%3D%20Expr.repeat(Expr.reshape(-Dt_r%2C%5B1%2CL%2Cm%5D)%2Cp%2C0)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Matrix%20mult(V%5ET%2CD)%2C%20stacked%20in%20m%2C%20dim%20%3D%20%5Bp%2CL%2Cm%5D.%0A%20%20%20%20%20%20%20%20%20%20%20%20Vt_D%20%3D%20Expr.stack(2%2C%5BVt%5B%3Ap%2C%3Aq%2Ci%5D.reshape(%5Bp%2Cq%5D)%20%40%20Matrix.sparse(D)%20for%20i%20in%20range(m)%5D)%0A%20%20%20%20%20%20%20%20%20%20%20%20h_23%20%3D%20Ax%20-%20Vt_D%20-%20Expr.repeat(Expr.reshape(Dt%40u%2C%5B1%2CL%2Cm%5D)%2Cp%2C0)%0A%20%20%20%20%20%20%20%20%20%20%20%20%23Second%20set%20of%20cones.%20Note%20the%20stacking%20in%204th%20dimension.%0A%20%20%20%20%20%20%20%20%20%20%20%20Cone2%20%3D%20M.constraint('Cone2'%2CExpr.stack(3%2C%5Bh_21%2Ch_22%2Ch_23%5D)%2CDomain.inPExpCone())%0A%20%20%20%20%20%20%20%20%20%20%20%20LinIneq2%20%3D%20M.constraint('LinIneq2'%2CExpr.sum(h_21%2C0)%20%2B%20Expr.neg(Dt_r)%20%3D%3D%200)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20%20%20%20%20return(M.primalObjValue()%2CM.getSolverDoubleInfo('optimizerTime'))%0A%20%20%20%20return%20(approx_2%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%23%23%20Approximation%20error%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20To%20evaluate%20the%20quality%20of%20the%20above%20approximations%2C%20Roos%20et.%20al.%20have%20used%20the%20following%20approximation%20error%20(percentage)%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20100%5Cbigg(%20%5Cfrac%7Be%5E%7Bc%5ET%20%5Chat%7Bx%7D%7D%7D%7Be%5E%7Bc%5ET%20x%5E*%7D%7D%20-1%5Cbigg)%0A%20%20%20%20%24%24%0A%0A%20%20%20%20We%20compare%20the%20objective%20values%20of%20the%20approximation%20to%20that%20of%20the%20exact%20solution.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20%23Making%20instances%20of%20the%20problem%20with%20n%3Dm%3D100%20and%20L%20ranging%20from%201%20to%20L_max.%0A%20%20%20%20%23Solve%20some%20initial%20ones%20exactly%20and%20others%20approximately%0A%20%20%20%20%23You%20can%20increase%20L_max%2C%20here%20we%20keep%20it%20small%20for%20the%20purpose%20of%20demonstration%0A%20%20%20%20L_max%20%3D%2010%0A%20%20%20%20t%2C%20tExact%2C%20n%2C%20m%20%3D%20L_max%2C%204%2C%20100%2C%20100%0A%20%20%20%20n_list%20%3D%20%5Bn%5D*t%0A%20%20%20%20m_list%20%3D%20%5Bm%5D*t%0A%20%20%20%20L_list%20%3D%20np.arange(1%2Ct%2B1%2C1)%0A%20%20%20%20return%20L_list%2C%20L_max%2C%20m_list%2C%20n_list%2C%20tExact%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22We%20will%20generate%20random%20matrices%2C%20with%20sparsity%20density%20%24s%20%3D%200.1%24.%20For%20every%20inequality%20constraint%20in%20the%20problem%2C%20%20there%20are%20%24p%20%3D%202%24%20such%20matrices.%20As%20stated%20before%2C%20%24m%3D100%24%20inequalities%20are%20involved.%20The%20dimension%20of%20the%20uncertain%20parameter%2C%20i.e.%20%24L%24%20(%24%5Czeta%20%5Cin%20%5CBbb%20R%5EL%24)%20ranges%20in%20%24%5B1%2CL_%7Bmax%7D%5D%24.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20Exact_sol%2C%0A%20%20%20%20Matrix%2C%0A%20%20%20%20SolutionError%2C%0A%20%20%20%20Sparse_make%2C%0A%20%20%20%20approx_1%2C%0A%20%20%20%20approx_2%2C%0A%20%20%20%20np%2C%0A%20%20%20%20tExact%2C%0A)%3A%0A%20%20%20%20def%20main(n_l%2Cm_l%2CL_l%2Cp%2Cs)%3A%0A%20%20%20%20%20%20%20%20Ex%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20App1%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20App2%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20for%20j%20in%20range(len(n_l))%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20print('%5Cn%20Iteration%3A%20%7B%7D'.format(j))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20Bt%20%3D%20%5B%5BMatrix.sparse(Sparse_make(n_l%5Bj%5D%2CL_l%5Bj%5D%2Cs%2Ci))%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Matrix.sparse(Sparse_make(n_l%5Bj%5D%2CL_l%5Bj%5D%2Cs%2Ci%2Bm_l%5Bj%5D))%5D%20for%20i%20in%20range(m_l%5Bj%5D)%5D%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Approximation%201%0A%20%20%20%20%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20print('%20%20%20Solve%20approximation%201')%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20App1.append(approx_1(Bt%2Cn_l%5Bj%5D%2Cm_l%5Bj%5D%2CL_l%5Bj%5D%2Cp))%0A%20%20%20%20%20%20%20%20%20%20%20%20except%20SolutionError%20as%20e%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20print('%20%20%20'%20%2B%20e.toString())%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Approximation%202%0A%20%20%20%20%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20print('%20%20%20Solve%20approximation%202')%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20App2.append(approx_2(Bt%2Cn_l%5Bj%5D%2Cm_l%5Bj%5D%2CL_l%5Bj%5D%2Cp))%0A%20%20%20%20%20%20%20%20%20%20%20%20except%20SolutionError%20as%20e%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20print('%20%20%20'%20%2B%20e.toString())%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20We%20only%20solve%20the%20exact%20problem%20for%20reasonably%20small%20j%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20j%3CtExact%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20try%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20print('%20%20%20Solve%20exact%20version')%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20Ex.append(Exact_sol(Bt%2Cn_l%5Bj%5D%2Cm_l%5Bj%5D%2CL_l%5Bj%5D%2Cp))%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20except%20SolutionError%20as%20e%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20print('%20%20%20'%20%2B%20e.toString())%0A%0A%20%20%20%20%20%20%20%20return(np.asarray(Ex)%2Cnp.asarray(App1)%2Cnp.asarray(App2))%0A%20%20%20%20return%20(main%2C)%0A%0A%0A%40app.cell%0Adef%20_(L_list%2C%20m_list%2C%20main%2C%20n_list)%3A%0A%20%20%20%20%23Sparsity%20density%20set%20to%200.1%20and%20number%20of%20terms%20in%20the%20log-sum-exp%20is%202%2C%20like%20in%20the%20paper.%0A%20%20%20%20s%20%3D%200.1%0A%20%20%20%20p%20%3D%202%0A%20%20%20%20Ex%2CApp1%2CApp2%20%3D%20main(n_list%2Cm_list%2CL_list%2Cp%2Cs)%0A%20%20%20%20return%20App1%2C%20App2%2C%20Ex%0A%0A%0A%40app.cell%0Adef%20_(App1%2C%20App2%2C%20Ex%2C%20L_max%2C%20np%2C%20plt%2C%20tExact)%3A%0A%20%20%20%20%23Plot%20for%20the%20computation%20time%20of%20the%20exact%20case.%0A%20%20%20%20L_exact_comp%20%3D%20np.arange(1%2Clen(Ex)%2B1%2C1)%0A%20%20%20%20L_comp%20%3D%20np.arange(1%2CL_max%2C1)%0A%20%20%20%20plt.plot(L_comp%2CApp1%5B0%3A(L_max-1)%2C1%5D%2Clabel%3D'Approx%201')%0A%20%20%20%20plt.plot(L_comp%2CApp2%5B0%3A(L_max-1)%2C1%5D%2Clabel%3D'Approx%202')%0A%20%20%20%20plt.plot(L_exact_comp%2CEx%5B%3A%2C1%5D%2Clabel%3D'Exact%20case')%0A%20%20%20%20plt.legend()%0A%20%20%20%20plt.xlabel('L')%0A%20%20%20%20plt.xticks(L_comp)%0A%20%20%20%20plt.ylabel('Optimizer%20Time%20(s)')%0A%20%20%20%20plt.title('Computation%20time%20for%20Exact%20Solution%20vs%20L')%0A%20%20%20%20plt.show()%0A%20%20%20%20plt.close()%0A%0A%20%20%20%20%23Plot%20for%20the%20computation%20times%20of%20the%20approximations.%0A%0A%20%20%20%20plt.plot(L_comp%2CApp1%5B0%3A(L_max-1)%2C1%5D%2Clabel%3D'Approx%201')%0A%20%20%20%20plt.plot(L_comp%2CApp2%5B0%3A(L_max-1)%2C1%5D%2Clabel%3D'Approx%202')%0A%20%20%20%20plt.legend()%0A%20%20%20%20plt.xlabel('L')%0A%20%20%20%20plt.xticks(L_comp)%0A%20%20%20%20plt.ylabel('Optimizer%20Time%20(s)')%0A%20%20%20%20plt.title('Computation%20time%20for%20approximations%20vs%20L')%0A%20%20%20%20plt.show()%0A%20%20%20%20plt.close()%0A%0A%20%20%20%20%23Error%20in%20the%20approximation%20vs%20L%0A%20%20%20%20L_err%20%3D%20np.arange(1%2CtExact%2B1%2C1)%0A%20%20%20%20error1%20%3D%20100*((np.exp(App1%5B0%3AtExact%2C0%5D)%2Fnp.exp(Ex%5B0%3AtExact%2C0%5D))-1)%0A%20%20%20%20error2%20%3D%20100*((np.exp(App2%5B0%3AtExact%2C0%5D)%2Fnp.exp(Ex%5B0%3AtExact%2C0%5D))-1)%0A%20%20%20%20plt.plot(L_err%2Cerror1%2Clabel%3D'Approx%201')%0A%20%20%20%20plt.plot(L_err%2Cerror2%2Clabel%3D'Approx%202')%0A%20%20%20%20plt.legend()%0A%20%20%20%20plt.xlabel('L')%0A%20%20%20%20plt.xticks(L_err)%0A%20%20%20%20plt.ylabel('Approximation%20error%20(%25)')%0A%20%20%20%20plt.title('Approximation%20error%20vs%20L')%0A%20%20%20%20plt.show()%0A%20%20%20%20plt.close()%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_()%3A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20return%20(mo%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20%23%20More%20results%0A%0A%20%20%20%20Results%20from%20an%20extended%20run%20of%20the%20above%20experiment%20are%20shown%20below.%20They%20were%20performed%20on%20a%20desktop%20with%2015.6%20GB%20of%20RAM%20and%20an%20Intel%24%5E%7B%5CcircledR%7D%24%20Core%24%5E%7B%5Ctext%7BTM%7D%7D%24%20i7-6770HQ%20CPU%20%40%202.6%20GHz%20%24%5Ctimes%24%208.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20_src%20%3D%20%22https%3A%2F%2Fraw.githubusercontent.com%2FMOSEK%2FTutorials%2Fmaster%2Fapprox-uncertain-ineq%2Ffullrun1.png%22%0A%0A%20%20%20%20mo.image(src%3D_src%2Crounded%3DTrue)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20_src2%20%3D%20%22https%3A%2F%2Fraw.githubusercontent.com%2FMOSEK%2FTutorials%2Fmaster%2Fapprox-uncertain-ineq%2Ffullrun2.png%22%0A%0A%20%20%20%20mo.image(src%3D_src2%2Crounded%3DTrue)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20_src3%20%3D%20%22https%3A%2F%2Fraw.githubusercontent.com%2FMOSEK%2FTutorials%2Fmaster%2Fapprox-uncertain-ineq%2Ffullrun3.png%22%0A%0A%20%20%20%20mo.image(src%3D_src3%2Crounded%3DTrue)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%3Ca%20rel%3D%22license%22%20href%3D%22http%3A%2F%2Fcreativecommons.org%2Flicenses%2Fby%2F4.0%2F%22%3E%3Cimg%20alt%3D%22Creative%20Commons%20License%22%20style%3D%22border-width%3A0%22%20src%3D%22https%3A%2F%2Fi.creativecommons.org%2Fl%2Fby%2F4.0%2F80x15.png%22%20%2F%3E%3C%2Fa%3E%3Cbr%20%2F%3EThis%20work%20is%20licensed%20under%20a%20%3Ca%20rel%3D%22license%22%20href%3D%22http%3A%2F%2Fcreativecommons.org%2Flicenses%2Fby%2F4.0%2F%22%3ECreative%20Commons%20Attribution%204.0%20International%20License%3C%2Fa%3E.%20The%20**MOSEK**%20logo%20and%20name%20are%20trademarks%20of%20%3Ca%20href%3D%22http%3A%2F%2Fmosek.com%22%3EMosek%20ApS%3C%2Fa%3E.%20The%20code%20is%20provided%20as-is.%20Compatibility%20with%20future%20release%20of%20**MOSEK**%20or%20the%20%60Fusion%20API%60%20are%20not%20guaranteed.%20For%20more%20information%20contact%20our%20%5Bsupport%5D(mailto%3Asupport%40mosek.com).%22%22%22)%0A%20%20%20%20return%0A%0A%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A
86f58cce28295ad1d69c82923c2cfb7d3d4108c0c49c8d323fe07bc1bf9127ea