control

import%20marimo%0A%0A__generated_with%20%3D%20%220.25.2%22%0Aapp%20%3D%20marimo.App(width%3D%22medium%22)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_()%3A%0A%20%20%20%20import%20sys%0A%20%20%20%20from%20pathlib%20import%20Path%0A%0A%20%20%20%20import%20marimo%20as%20mo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20scipy.sparse%20as%20spa%0A%0A%20%20%20%20start%20%3D%20Path(__file__).resolve().parent%20if%20%22__file__%22%20in%20globals()%20else%20Path.cwd().resolve()%0A%20%20%20%20for%20candidate%20in%20(start%2C%20*start.parents)%3A%0A%20%20%20%20%20%20%20%20if%20(candidate%20%2F%20%22fomtools%22).is_dir()%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20repo_root%20%3D%20candidate%0A%20%20%20%20%20%20%20%20%20%20%20%20break%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20repo_root%20%3D%20start%0A%0A%20%20%20%20if%20str(repo_root)%20not%20in%20sys.path%3A%0A%20%20%20%20%20%20%20%20sys.path.insert(0%2C%20str(repo_root))%0A%0A%20%20%20%20from%20docs.tutorials.helper_functions.control%20import%20(%0A%20%20%20%20%20%20%20%20admm_step%2C%0A%20%20%20%20%20%20%20%20build_admm_lhs%2C%0A%20%20%20%20%20%20%20%20build_control_qp%2C%0A%20%20%20%20%20%20%20%20estimate_control_pepit_radius%2C%0A%20%20%20%20%20%20%20%20make_control_dynamics%2C%0A%20%20%20%20)%0A%20%20%20%20from%20fomtools.proxable_functions%20import%20IndBox%0A%0A%20%20%20%20np.set_printoptions(precision%3D5)%20%20%23%20Print%20few%20decimal%20places%0A%20%20%20%20np.set_printoptions(suppress%3DTrue)%20%20%23%20Suppress%20scientific%20notation%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20IndBox%2C%0A%20%20%20%20%20%20%20%20admm_step%2C%0A%20%20%20%20%20%20%20%20build_admm_lhs%2C%0A%20%20%20%20%20%20%20%20build_control_qp%2C%0A%20%20%20%20%20%20%20%20estimate_control_pepit_radius%2C%0A%20%20%20%20%20%20%20%20make_control_dynamics%2C%0A%20%20%20%20%20%20%20%20mo%2C%0A%20%20%20%20%20%20%20%20np%2C%0A%20%20%20%20%20%20%20%20plt%2C%0A%20%20%20%20%20%20%20%20spa%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%20Problem%20setup%3A%20Constrained%20optimal%20control%0A%0A%20%20%20%20%5BBack%20to%20the%20verification%20tutorials%5D(.%2F)%0A%0A%20%20%20%20In%20this%20tutorial%2C%20we%20verify%20an%20OSQP-style%20ADMM%20method%20for%20a%20finite-horizon%20constrained%20optimal%20control%20problem.%20The%20system%20has%20state%20%24s_t%20%5Cin%20%5Cmathbf%7BR%7D%5E%7Bn_s%7D%24%2C%20input%20%24u_t%20%5Cin%20%5Cmathbf%7BR%7D%5E%7Bn_u%7D%24%2C%20and%20initial%20state%20parameter%20%24x%24.%0A%0A%20%20%20%20%5Cbegin%7Barray%7D%7Blll%7D%0A%20%20%20%20%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20s_T%5ET%20Q_T%20s_T%20%2B%20%5Csum_%7Bt%3D1%7D%5E%7BT-1%7D%20s_t%5ET%20Q%20s_t%20%2B%20u_t%5ET%20R%20u_t%20%5C%5C%0A%20%20%20%20%20%20%20%20%5Ctext%7Bsubject%20to%7D%20%26%20s_%7Bt%2B1%7D%20%3D%20A%20s_t%20%2B%20B%20u_t%2C%20%26%20t%3D1%2C%5Cdots%2CT-1%20%5C%5C%0A%20%20%20%20%20%20%20%20%26%20s_%7B%5Cmin%7D%20%5Cle%20s_t%20%5Cle%20s_%7B%5Cmax%7D%2C%20%26%20t%3D1%2C%5Cdots%2CT%20%5C%5C%0A%20%20%20%20%20%20%20%20%26%20u_%7B%5Cmin%7D%20%5Cle%20u_t%20%5Cle%20u_%7B%5Cmax%7D%2C%20%26%20t%3D1%2C%5Cdots%2CT-1%20%5C%5C%0A%20%20%20%20%20%20%20%20%26%20s_1%20%3D%20x%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20with%20optimization%20variable%20%24z%20%3D%20(u_1%2C%20%5Cdots%2C%20u_%7BT-1%7D)%20%5Cin%20%5Cmathbf%7BR%7D%5E%7B(T-1)%20n_u%7D%24%20and%20parameter%20%24x%20%5Cin%20%5Cmathbf%7BR%7D%5E%7Bn_s%7D%24.%0A%20%20%20%20The%20verification%20problem%20searches%20over%20a%20bounded%20set%20of%20initial%20ADMM%20iterates%20and%2C%20depending%20on%20the%20experiment%2C%20either%20fixes%20%24x%24%20or%20lets%20it%20vary%20over%20the%20prescribed%20box.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(make_control_dynamics%2C%20np%2C%20spa)%3A%0A%20%20%20%20%23%20problem%20data%3A%20state%20dimension%20n_s%20%3D%2012%2C%20input%20dimension%20n_u%20%3D%204%0A%20%20%20%20A_dym%2C%20B_dym%20%3D%20make_control_dynamics()%0A%20%20%20%20(n_s%2C%20n_u)%20%3D%20B_dym.shape%0A%0A%20%20%20%20T%20%3D%203%20%20%20%20%20%20%20%20%23%20finite%20horizon%20t%3D1%2C2%2C...%2CT%0A%0A%20%20%20%20Q%20%3D%20spa.diags(%5B0.%2C%200.%2C%2010.%2C%2010.%2C%2010.%2C%2010.%2C%200.%2C%200.%2C%200.%2C%205.%2C%205.%2C%205.%5D)%20%23%20state%20cost%20matrix%20Q_t%20%3D%20Q_T%20%3D%20Q%0A%20%20%20%20R%20%3D%200.1*spa.eye(n_u)%20%23%20input%20cost%20matrix%0A%0A%20%20%20%20%23%20State%20constraints%20are%20unconstrained%20in%20this%20tutorial.%0A%0A%20%20%20%20%23%20initial%20state%20s_init%20as%20parameter%20x%3A%20the%20set%20X%20%3D%20%5Bx_lb%2C%20x_ub%5D%0A%20%20%20%20x_lb%20%3D%20np.array(%5B-np.pi%2F6%2C%20-np.pi%2F6%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%5D)%0A%20%20%20%20x_ub%20%3D%20np.array(%5B%20np.pi%2F6%2C%20%20np.pi%2F6%2C%201.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%5D)%0A%20%20%20%20x_init%20%3D%202.0*np.array(%5B1.%2C%201.%2C%201.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%2C%200.%5D)%20%23%20initial%20x%20when%20we%20fix%20it%0A%0A%20%20%20%20%23%20input%20constraint%0A%20%20%20%20u_min%20%3D%20-0.99*np.ones(n_u)%0A%20%20%20%20u_max%20%3D%20%202.41*np.ones(n_u)%0A%20%20%20%20u_stack_min%20%3D%20np.tile(u_min%2C%20T%20-%201)%0A%20%20%20%20u_stack_max%20%3D%20np.tile(u_max%2C%20T%20-%201)%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20A_dym%2C%0A%20%20%20%20%20%20%20%20B_dym%2C%0A%20%20%20%20%20%20%20%20Q%2C%0A%20%20%20%20%20%20%20%20R%2C%0A%20%20%20%20%20%20%20%20T%2C%0A%20%20%20%20%20%20%20%20n_s%2C%0A%20%20%20%20%20%20%20%20n_u%2C%0A%20%20%20%20%20%20%20%20u_stack_max%2C%0A%20%20%20%20%20%20%20%20u_stack_min%2C%0A%20%20%20%20%20%20%20%20x_init%2C%0A%20%20%20%20%20%20%20%20x_lb%2C%0A%20%20%20%20%20%20%20%20x_ub%2C%0A%20%20%20%20)%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20%23%20ADMM%20parameters%0A%20%20%20%20sigma%20%3D%201e-4%0A%20%20%20%20rho%20%3D%201.%20%20%23%200.1%2C%201%2C%2010%0A%20%20%20%20return%20rho%2C%20sigma%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Quadratic%20program%20form%0A%0A%20%20%20%20Now%2C%20consider%20the%20equivalent%20problem%20of%20the%20form%3A%0A%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20(1%2F2)%20z%5ET%20P%20z%20%2B%20q(x)%5ET%20z%20%5C%5C%0A%20%20%20%20%20%20%20%20%5Ctext%7Bsubject%20to%7D%20%26%20l(x)%20%5Cle%20Mz%20%5Cle%20u(x)%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20where%20%24z%20%3D%20(u_1%2C%20%5Cdots%2C%20u_%7BT-1%7D)%24.%20Note%20that%20for%20%24t%3D0%2C%5Cdots%2CT-1%24%2C%0A%0A%20%20%20%20%24%24%20s_%7Bt%2B1%7D%20%3D%20A%20s_t%20%2B%20B%20u_t%20%3D%20A%5Et%20x%20%2B%20%5Csum_%7Bm%3D1%7D%5E%7Bt%7D%20A%5E%7Bt-m%7D%20B%20u_m%20%3D%20A%5Et%20x%20%2B%20%5Cbegin%7Bbmatrix%7D%20A%5E%7Bt-1%7D%20B%20%26%20%5Cdots%20%26%20B%20%26%200%20%26%20%5Cdots%20%26%200%20%5Cend%7Bbmatrix%7D%20z.%20%24%24%0A%0A%20%20%20%20In%20verification%20problem%20(VP)%20implementation%2C%20we%20use%20the%20optimization%20variable%20%60xz_stack%60%20defined%20as%20%24(x%2C%20u_1%2C%20%5Cdots%2C%20u_%7BT-1%7D)%20%5Cin%20%5Cmathbf%7BR%7D%5E%7Bn_s%20%2B%20(T-1)%20n_u%7D%24.%0A%20%20%20%20This%20makes%20the%20initial%20state%20constraints%20and%20input%20box%20constraints%20part%20of%20one%20stacked%20box%20constraint.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(A_dym%2C%20B_dym%2C%20Q%2C%20R%2C%20T%2C%20build_control_qp%2C%20np)%3A%0A%20%20%20%20M%2C%20P%2C%20q%2C%20n_xz%20%3D%20build_control_qp(A_dym%2C%20B_dym%2C%20Q%2C%20R%2C%20T)%0A%20%20%20%20z0_lb%2C%20z0_ub%20%3D%20-np.ones(n_xz)%2C%20np.ones(n_xz)%0A%20%20%20%20v0_lb%2C%20v0_ub%20%3D%20np.zeros(n_xz)%2C%20np.ones(n_xz)%0A%20%20%20%20return%20M%2C%20P%2C%20n_xz%2C%20q%2C%20v0_lb%2C%20v0_ub%2C%20z0_lb%2C%20z0_ub%0A%0A%0A%40app.cell%0Adef%20_(M%2C%20build_admm_lhs%2C%20rho%2C%20sigma)%3A%0A%20%20%20%20def%20build_lhs(P)%3A%0A%20%20%20%20%20%20%20%20%22%22%22OSQP%20linear-solve%20matrix%20%60%60P%20%2B%20sigma%20I%20%2B%20rho%20M%5ET%20M%60%60%20(sparse).%22%22%22%0A%20%20%20%20%20%20%20%20return%20build_admm_lhs(M%2C%20P%2C%20rho%2C%20sigma)%0A%0A%20%20%20%20return%20(build_lhs%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20First-order%20method%3A%20OSQP-style%20ADMM%0A%0A%20%20%20%20We%20solve%20the%20boxed%20quadratic%20program%20with%20the%20splitting%20used%20by%20OSQP%3A%0A%0A%20%20%20%20%5Cbegin%7Barray%7D%7Brcl%7D%0A%20%20%20%20%20%20%20%20w%5E%7Bk%2B1%7D%20%26%3D%26%20%5CPi_%7B%5Bl(x)%2C%20u(x)%5D%7D(v%5Ek)%2C%20%5C%5C%0A%20%20%20%20%20%20%20%20z%5E%7Bk%2B1%7D%20%26%3D%26%20(P%20%2B%20%5Csigma%20I%20%2B%20%5Crho%20M%5ET%20M)%5E%7B-1%7D%0A%20%20%20%20%20%20%20%20%20%20%20%20%5Cleft(%5Csigma%20z%5Ek%20-%20q%20%2B%20%5Crho%20M%5ET(2w%5E%7Bk%2B1%7D%20-%20v%5Ek)%5Cright)%2C%20%5C%5C%0A%20%20%20%20%20%20%20%20v%5E%7Bk%2B1%7D%20%26%3D%26%20v%5Ek%20%2B%20Mz%5E%7Bk%2B1%7D%20-%20w%5E%7Bk%2B1%7D.%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20To%20compare%20all%20verification%20methods%20in%20the%20same%20units%2C%20we%20measure%20the%0A%20%20%20%20squared%20Euclidean%20fixed-point%20residual%20of%20the%20OSQP%20state%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%20%20%20%20%5Cleft%5C%7C%20z%5E%7Bk%2B1%7D%20-%20z%5Ek%20%5Cright%5C%7C_2%5E2%20%2B%20%5Cleft%5C%7C%20v%5E%7Bk%2B1%7D%20-%20v%5Ek%20%5Cright%5C%7C_2%5E2.%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20projection%20%24%5CPi_%7B%5Bl(x)%2C%20u(x)%5D%7D%24%20forces%20the%20state%20coordinates%20to%20equal%0A%20%20%20%20%24x%24%20and%20clips%20the%20input%20coordinates%20to%20%24%5Bu_%7B%5Cmin%7D%2C%20u_%7B%5Cmax%7D%5D%24.%20Because%20the%0A%20%20%20%20simulator%20cannot%20project%20onto%20the%20resulting%20zero-width%20state%20box%2C%20each%20mode%0A%20%20%20%20passes%20the%20appropriate%20projection%20into%20the%20shared%20%60admm_step%60%3A%20the%0A%20%20%20%20simulator%20clips%20the%20input%20block%20with%20%60satlin%60%20and%20stacks%20%24x%24%20on%20top%2C%20while%0A%20%20%20%20the%20MIP%20verification%20backend%20projects%20onto%20an%20%60IndBox%60%20whose%0A%20%20%20%20bounds%20pin%20the%20state%20block%20to%20%24x%24.%0A%20%20%20%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%0A%20%20%20%20%23%20Verification%20Problem%0A%0A%20%20%20%20We%20now%20describe%20how%20to%20use%20%60fomtools%60%20to%20verify%20the%20worst-case%20convergence%20rate%20of%20the%20OSQP-style%20ADMM%20method%20for%20constrained%20optimal%20control.%20We%20compare%20randomly%20sampled%20simulations%2C%20exact%20MIP%20verification%2C%20and%20a%20PEPit%20bound%2C%20up%20to%20%24K%24%20iterations.%20The%20initial%20state%20parameter%20%24x%24%20can%20either%20be%20fixed%20at%20%24x_%5Ctext%7Binit%7D%24%20or%20allowed%20to%20vary%20over%20its%20prescribed%20box.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_()%3A%0A%20%20%20%20K%20%3D%2020%20%20%23%20Maximum%20iteration%20number%0A%20%20%20%20fix_initial_x%20%3D%20True%20%20%23%20False%20lets%20x%20vary%20over%20%5Bx_lb%2C%20x_ub%5D.%0A%20%20%20%20return%20K%2C%20fix_initial_x%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Simulator%20class%0A%0A%20%20%20%20We%20first%20simulate%20the%20ADMM%20algorithm%20by%20sampling%20the%20initial%20iterates%0A%20%20%20%20%24(z%5E0%2Cv%5E0)%24%20over%20their%20prescribed%20boxes.%20The%20parameter%20%24x%24%20is%20either%20fixed%0A%20%20%20%20at%20%24x_%5Ctext%7Binit%7D%24%20or%20sampled%20from%20%24%5Bx_%5Cmathrm%7Blb%7D%2Cx_%5Cmathrm%7Bub%7D%5D%24%2C%0A%20%20%20%20according%20to%20%60fix_initial_x%60%3B%20the%20MIP%20and%20radius%20computation%20use%20the%20same%0A%20%20%20%20choice.%20The%20samples%20are%20drawn%20once%20below%20and%20reused%20by%20every%20simulation%0A%20%20%20%20cell.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20fix_initial_x%2C%0A%20%20%20%20n_s%2C%0A%20%20%20%20n_xz%2C%0A%20%20%20%20np%2C%0A%20%20%20%20v0_lb%2C%0A%20%20%20%20v0_ub%2C%0A%20%20%20%20x_init%2C%0A%20%20%20%20x_lb%2C%0A%20%20%20%20x_ub%2C%0A%20%20%20%20z0_lb%2C%0A%20%20%20%20z0_ub%2C%0A)%3A%0A%20%20%20%20n_samples%20%3D%2010000%0A%20%20%20%20_rng%20%3D%20np.random.default_rng(0)%0A%0A%20%20%20%20if%20fix_initial_x%3A%0A%20%20%20%20%20%20%20%20x_samples%20%3D%20np.repeat(x_init%5Bnp.newaxis%2C%20%3A%5D%2C%20n_samples%2C%20axis%3D0)%0A%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20x_samples%20%3D%20_rng.uniform(low%3Dx_lb%2C%20high%3Dx_ub%2C%20size%3D(n_samples%2C%20n_s))%0A%20%20%20%20z0_samples%20%3D%20_rng.uniform(low%3Dz0_lb%2C%20high%3Dz0_ub%2C%20size%3D(n_samples%2C%20n_xz))%0A%20%20%20%20v0_samples%20%3D%20_rng.uniform(low%3Dv0_lb%2C%20high%3Dv0_ub%2C%20size%3D(n_samples%2C%20n_xz))%0A%20%20%20%20return%20v0_samples%2C%20x_samples%2C%20z0_samples%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20The%20simulation%20is%20run%20through%20the%20following%20function.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20M%2C%0A%20%20%20%20admm_step%2C%0A%20%20%20%20build_lhs%2C%0A%20%20%20%20n_s%2C%0A%20%20%20%20n_xz%2C%0A%20%20%20%20q%2C%0A%20%20%20%20rho%2C%0A%20%20%20%20sigma%2C%0A%20%20%20%20u_stack_max%2C%0A%20%20%20%20u_stack_min%2C%0A)%3A%0A%20%20%20%20import%20fomtools%20as%20ft%0A%0A%20%20%20%20def%20maximize_residual(vp%2C%20residual%2C%20ord)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Maximize%20%60%60%7C%7Cresidual%7C%7C_ord%60%60%2C%20squared%20for%20%60%60ord%3D2%60%60%20(%60%60ft.sum_squares%60%60).%22%22%22%0A%20%20%20%20%20%20%20%20vp.maximize(ft.sum_squares(residual)%20if%20ord%20%3D%3D%202%20else%20ft.norm(residual%2C%20ord))%0A%0A%20%20%20%20def%20admm_simulation(x_samples%2C%20z0_samples%2C%20v0_samples%2C%20P%2C%20K%2C%20ord%3D2)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Simulate%20ADMM%20on%20sampled%20instances.%0A%0A%20%20%20%20%20%20%20%20Builds%20a%20simulator-mode%20%3Aclass%3A%60Verifier%60%2C%20registers%20the%20initial-state%0A%20%20%20%20%20%20%20%20parameter%20%60%60x%60%60%20and%20the%20initial%20iterates%20%60%60(z0%2C%20v0)%60%60%20as%20sampled%0A%20%20%20%20%20%20%20%20parameters%2C%20rolls%20out%20the%20recurrence%20with%20%60%60admm_step%60%60%2C%20and%20simulates%0A%20%20%20%20%20%20%20%20via%20%60%60vp.solve()%60%60.%0A%20%20%20%20%20%20%20%20The%20projection%20clips%20the%20input%20block%20with%20%60%60satlin%60%60%20and%20stacks%20%60%60x%60%60%0A%20%20%20%20%20%20%20%20on%20top.%20Each%20%60%60*_samples%60%60%20array%20has%20shape%20%60%60(n_samples%2C%20dim)%60%60%20(a%0A%20%20%20%20%20%20%20%20single-row%20array%20replays%20one%20instance)%2C%20and%20the%20returned%20%60%60residuals%60%60%0A%20%20%20%20%20%20%20%20has%20shape%20%60%60(n_samples%2C%20K)%60%60.%20Summary%20statistics%20(%60%60mean%60%60%2C%20%60%60min%60%60%2C%0A%20%20%20%20%20%20%20%20%60%60max%60%60%2C%20...)%20are%20available%20on%20the%20returned%20instance%20%60%60vp%60%60%20via%0A%20%20%20%20%20%20%20%20%60%60vp%5B...%5D%60%60.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20vp%20%3D%20ft.Verifier(method%3D'simulator'%2C%20obbt_solver%3DNone)%0A%20%20%20%20%20%20%20%20x_param%20%3D%20vp.parameter(ndim%3Dn_s%2C%20samples%3Dx_samples)%0A%20%20%20%20%20%20%20%20z%20%3D%20vp.parameter(n_xz%2C%20samples%3Dz0_samples)%0A%20%20%20%20%20%20%20%20v%20%3D%20vp.parameter(n_xz%2C%20samples%3Dv0_samples)%0A%0A%20%20%20%20%20%20%20%20def%20project(vp%2C%20v_prev)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20w_control%20%3D%20vp.satlin(v_prev%5Bn_s%3A%5D%2C%20l%3Du_stack_min%2C%20u%3Du_stack_max)%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20ft.stack(%5Bx_param%2C%20w_control%5D)%0A%0A%20%20%20%20%20%20%20%20lhs_mat%20%3D%20build_lhs(P).toarray()%0A%20%20%20%20%20%20%20%20for%20_%20in%20range(K)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20z%2C%20v%2C%20residual%20%3D%20admm_step(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20vp%2C%20M%2C%20q%2C%20rho%2C%20sigma%2C%20z%2C%20v%2C%20project%2C%20lhs_mat%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20maximize_residual(vp%2C%20residual%2C%20ord)%0A%0A%20%20%20%20%20%20%20%20residuals%20%3D%20vp.solve()%0A%20%20%20%20%20%20%20%20return%20vp%2C%20residuals%0A%0A%20%20%20%20return%20admm_simulation%2C%20ft%2C%20maximize_residual%0A%0A%0A%40app.cell%0Adef%20_(K%2C%20P%2C%20admm_simulation%2C%20np%2C%20v0_samples%2C%20x_samples%2C%20z0_samples)%3A%0A%20%20%20%20_vp%2C%20_residuals%20%3D%20admm_simulation(x_samples%2C%20z0_samples%2C%20v0_samples%2C%20P%2C%20K%2C%20ord%3D2)%0A%0A%20%20%20%20print(%0A%20%20%20%20%20%20%20%20f%22Summary%20statistics%3A%5Cn%22%0A%20%20%20%20%20%20%20%20f%22mean%3A%20%7B_vp%5B'mean'%5D%3A.6e%7D%5Cn%22%0A%20%20%20%20%20%20%20%20f%22min%3A%20%7B_vp%5B'min'%5D%3A.6e%7D%5Cn%22%0A%20%20%20%20%20%20%20%20f%22max%3A%20%7B_vp%5B'max'%5D%3A.6e%7D%5Cn%22%0A%20%20%20%20%20%20%20%20f%22median%3A%20%7B_vp%5B'median'%5D%3A.6e%7D%5Cn%22%0A%20%20%20%20)%0A%0A%20%20%20%20simulation_residuals_admm%20%3D%20np.mean(_residuals%2C%20axis%3D0)%0A%20%20%20%20%23%20Quantile%201.0%20gives%20the%20in-sample%20worst%20case%3B%20indices%20are%20%5B90%25%2C%2050%25%2C%2010%25%2C%20max%5D.%0A%20%20%20%20simulation_residual_quantiles_admm%20%3D%20np.quantile(_residuals%2C%20%5B0.90%2C%200.50%2C%200.10%2C%201.0%5D%2C%20axis%3D0)%0A%20%20%20%20return%20simulation_residual_quantiles_admm%2C%20simulation_residuals_admm%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Algorithm%20verification%20through%20verification%20problem%20(VP)%0A%0A%20%20%20%20%23%23%23%20MIP-VP%0A%20%20%20%20Maximizing%20a%20squared%20Euclidean%20norm%20would%20make%20the%20verification%20problem%20a%0A%20%20%20%20nonconvex%20MIQP.%20Instead%2C%20the%20MIP%20maximizes%20the%20infinity%20norm%20of%20the%20same%0A%20%20%20%20OSQP-state%20residual%0A%20%20%20%20%24%24%0A%20%20%20%20%20%20%20%20r%5Ek%3D(z%5E%7Bk%2B1%7D-z%5Ek%2C%5C%3Bv%5E%7Bk%2B1%7D-v%5Ek).%0A%20%20%20%20%24%24%0A%20%20%20%20Since%20%24r%5Ek%5Cin%5Cmathbf%7BR%7D%5E%7B2n_%7Bxz%7D%7D%24%2C%20the%20norm%20inequality%0A%20%20%20%20%24%24%0A%20%20%20%20%20%20%20%20%5C%7Cr%5Ek%5C%7C_2%5E2%20%5Cle%202n_%7Bxz%7D%5C%7Cr%5Ek%5C%7C_%5Cinfty%5E2%0A%20%20%20%20%24%24%0A%20%20%20%20converts%20the%20exact%20MIP%20value%20into%20a%20valid%20upper%20bound%20in%20the%20squared%0A%20%20%20%20Euclidean%20units%20used%20by%20PEP%20and%20the%20simulations.%20The%20verification%20below%0A%20%20%20%20rolls%20out%20the%20same%20%60admm_step%60%2C%20projecting%20onto%20an%20%60IndBox%60%20whose%20bounds%0A%20%20%20%20pin%20the%20state%20block%20to%20%24x%24.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20IndBox%2C%0A%20%20%20%20M%2C%0A%20%20%20%20T%2C%0A%20%20%20%20admm_step%2C%0A%20%20%20%20build_lhs%2C%0A%20%20%20%20fix_initial_x%2C%0A%20%20%20%20ft%2C%0A%20%20%20%20maximize_residual%2C%0A%20%20%20%20n_s%2C%0A%20%20%20%20n_u%2C%0A%20%20%20%20n_xz%2C%0A%20%20%20%20np%2C%0A%20%20%20%20q%2C%0A%20%20%20%20rho%2C%0A%20%20%20%20sigma%2C%0A%20%20%20%20u_stack_max%2C%0A%20%20%20%20u_stack_min%2C%0A%20%20%20%20v0_lb%2C%0A%20%20%20%20v0_ub%2C%0A%20%20%20%20x_init%2C%0A%20%20%20%20x_lb%2C%0A%20%20%20%20x_ub%2C%0A%20%20%20%20z0_lb%2C%0A%20%20%20%20z0_ub%2C%0A)%3A%0A%20%20%20%20def%20admm_verification(K%2C%20method%2C%20P%2C%20ord%3Dnp.inf%2C%20**verifier_kwargs)%3A%0A%20%20%20%20%20%20%20%20%22%22%22Solve%20the%20optimal-control%20ADMM%20verification%20problem%20for%20one%20backend.%0A%0A%20%20%20%20%20%20%20%20Builds%20a%20%60%60Verifier%60%60%20in%20%60%60method%60%60%20mode%2C%20registers%20the%20initial-state%0A%20%20%20%20%20%20%20%20parameter%20(fixed%20at%20%60%60x_init%60%60%20or%20searched%20over%20%60%60%5Bx_lb%2C%20x_ub%5D%60%60%20per%0A%20%20%20%20%20%20%20%20%60%60fix_initial_x%60%60)%20and%20the%20initial%20iterate%20boxes%2C%20and%20rolls%20out%0A%20%20%20%20%20%20%20%20%60%60admm_step%60%60%20projecting%20onto%20the%20stacked%20%60%60IndBox%60%60.%20Returns%0A%20%20%20%20%20%20%20%20%60%60(vp%2C%20residuals)%60%60%20with%20one%20worst-case%20residual%20per%20iteration%3B%20the%20MIP%0A%20%20%20%20%20%20%20%20uses%20the%20infinity-norm%20objective%2C%20so%20the%20caller%20converts%20it%20to%20squared%0A%20%20%20%20%20%20%20%20Euclidean%20units.%0A%20%20%20%20%20%20%20%20%22%22%22%0A%20%20%20%20%20%20%20%20vp%20%3D%20ft.Verifier(method%3Dmethod%2C%20**verifier_kwargs)%0A%0A%20%20%20%20%20%20%20%20if%20fix_initial_x%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20x_param%20%3D%20vp.parameter(n_s%2C%20lb%3Dx_init%2C%20ub%3Dx_init)%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20x_param%20%3D%20vp.parameter(n_s%2C%20lb%3Dx_lb%2C%20ub%3Dx_ub)%0A%0A%20%20%20%20%20%20%20%20%23%20l(x)%20%3C%3D%20M%40z%20%3C%3D%20u(x)%3A%20pin%20the%20state%20block%20to%20x%2C%20clip%20the%20input%20block.%0A%20%20%20%20%20%20%20%20l_param%20%3D%20vp.parameter((T-1)*n_u%2C%20lb%3Du_stack_min%2C%20ub%3Du_stack_min)%0A%20%20%20%20%20%20%20%20u_param%20%3D%20vp.parameter((T-1)*n_u%2C%20lb%3Du_stack_max%2C%20ub%3Du_stack_max)%0A%20%20%20%20%20%20%20%20l_x%20%3D%20ft.stack(%5Bx_param%2C%20l_param%5D)%0A%20%20%20%20%20%20%20%20u_x%20%3D%20ft.stack(%5Bx_param%2C%20u_param%5D)%0A%20%20%20%20%20%20%20%20indbox%20%3D%20IndBox(l%3Dl_x%2C%20u%3Du_x)%0A%0A%20%20%20%20%20%20%20%20def%20project(vp%2C%20v_prev)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20vp.prox(indbox%2C%20v_prev)%0A%0A%20%20%20%20%20%20%20%20z%20%3D%20vp.iterate_trajectory(K%20%2B%201%2C%20collapse_expr%3D'auto')%0A%20%20%20%20%20%20%20%20v%20%3D%20vp.iterate_trajectory(K%20%2B%201%2C%20collapse_expr%3D'auto')%0A%20%20%20%20%20%20%20%20z%5B0%5D%20%3D%20vp.iterate(n_xz%2C%20lb%3Dz0_lb%2C%20ub%3Dz0_ub)%0A%20%20%20%20%20%20%20%20v%5B0%5D%20%3D%20vp.iterate(n_xz%2C%20lb%3Dv0_lb%2C%20ub%3Dv0_ub)%0A%20%20%20%20%20%20%20%20lhs_mat%20%3D%20build_lhs(P)%0A%0A%20%20%20%20%20%20%20%20residuals%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20for%20k%20in%20range(1%2C%20K%20%2B%201)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20z%2C%20v%2C%20residual%20%3D%20admm_step(%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20vp%2C%20M%2C%20q%2C%20rho%2C%20sigma%2C%20z%2C%20v%2C%20project%2C%20lhs_mat%2C%20k%3Dk%0A%20%20%20%20%20%20%20%20%20%20%20%20)%0A%20%20%20%20%20%20%20%20%20%20%20%20maximize_residual(vp%2C%20residual%2C%20ord)%0A%20%20%20%20%20%20%20%20%20%20%20%20residuals.append(vp.solve())%0A%20%20%20%20%20%20%20%20return%20vp%2C%20residuals%0A%0A%20%20%20%20return%20(admm_verification%2C)%0A%0A%0A%40app.cell%0Adef%20_(K%2C%20P%2C%20admm_verification%2C%20n_xz%2C%20np)%3A%0A%20%20%20%20%23%20Solve%20VP%20for%20optimal%20control%20example%20(infinity-norm%20objective).%0A%20%20%20%20_%2C%20all_residuals_inf%20%3D%20admm_verification(%0A%20%20%20%20%20%20%20%20K%2C%20'mip'%2C%20P%2C%20ord%3Dnp.inf%2C%20solver%3D'gurobi'%2C%0A%20%20%20%20)%0A%0A%20%20%20%20%23%20Convert%20the%20infinity-norm%20guarantee%20to%20squared%20Euclidean%20units.%0A%20%20%20%20all_residuals%20%3D%202%20*%20n_xz%20*%20np.square(all_residuals_inf)%0A%20%20%20%20print(f%22infinity-norm%20MIP%20bounds%3A%20%7Ball_residuals_inf%7D%22)%0A%20%20%20%20print(f%22squared-Euclidean%20upper%20bounds%3A%20%7Ball_residuals%7D%22)%0A%20%20%20%20return%20(all_residuals%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%23%20PEP-obtained%20worst-case%20convergence%20rate%0A%0A%20%20%20%20We%20compare%20the%20MIP%20result%20with%20a%20PEP%20that%20uses%20the%20same%20OSQP%20recurrence%2C%0A%20%20%20%20including%20the%20projection-first%20ordering%20and%20the%20parameters%20%24%5Csigma%24%20and%0A%20%20%20%20%24%5Crho%24.%20For%20%24M%3DI%24%2C%20the%20linear%20solve%20is%20the%20proximal%20step%0A%20%20%20%20%24%24%0A%20%20%20%20%20%20%20%20z%5E%7Bk%2B1%7D%0A%20%20%20%20%20%20%20%20%3D%20%5Cmathbf%7Bprox%7D_%7Bf%2F(%5Csigma%2B%5Crho)%7D%0A%20%20%20%20%20%20%20%20%5Cleft(%0A%20%20%20%20%20%20%20%20%20%20%20%20(%5Csigma%20z%5Ek%2B%5Crho(2w%5E%7Bk%2B1%7D-v%5Ek))%2F(%5Csigma%2B%5Crho)%0A%20%20%20%20%20%20%20%20%5Cright).%0A%20%20%20%20%24%24%0A%0A%20%20%20%20PEP%20still%20considers%20the%20larger%20class%20of%20all%20%24L%24-smooth%2C%0A%20%20%20%20%24%5Cmu%24-strongly-convex%20functions%20and%20all%20convex%20indicator%20sets%3A%0A%0A%20%20%20%20%5Cbegin%7Barray%7D%7Bll%7D%0A%20%20%20%20%20%20%20%20%5Ctext%7Bminimize%7D%20%26%20f(z)%20%2B%20g(z)%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%0A%20%20%20%20so%20it%20does%20not%20retain%20the%20exact%20matrix%20%24P%24%2C%20dimension%2C%20or%20box%20geometry.%20We%0A%20%20%20%20show%20both%20the%20bound%20conditioned%20on%20the%20empirical%20radius%20above%20and%20a%0A%20%20%20%20rate-shape%20curve%20calibrated%20to%20match%20the%20norm-converted%20MIP%20upper%20bound%20at%0A%20%20%20%20%24K%3D1%24.%20The%20calibrated%20curve%20is%20not%20an%20independent%20absolute%20guarantee.%0A%0A%20%20%20%20For%20the%20radius-conditioned%20PEP%20comparison%2C%20we%20estimate%20the%20OSQP-state%0A%20%20%20%20initial%20radius%0A%20%20%20%20%24%24%0A%20%20%20%20%20%20%20%20%5C%7Cz%5E0-z%5E%5Cstar%5C%7C_2%5E2%2B%5C%7Cv%5E0-v%5E%5Cstar%5C%7C_2%5E2%20%5Cle%20R%5E2%0A%20%20%20%20%24%24%0A%20%20%20%20from%20the%20same%20initial-state%20parameter%20set%20used%20by%20the%20MIP.%20Because%20the%0A%20%20%20%20maximum%20below%20is%20taken%20over%20random%20initial%20iterates%2C%20it%20is%20an%20empirical%0A%20%20%20%20radius%20rather%20than%20a%20certificate%20for%20every%20point%20in%20the%20initialization%0A%20%20%20%20boxes.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(%0A%20%20%20%20K%2C%0A%20%20%20%20M%2C%0A%20%20%20%20P%2C%0A%20%20%20%20all_residuals%2C%0A%20%20%20%20build_lhs%2C%0A%20%20%20%20estimate_control_pepit_radius%2C%0A%20%20%20%20np%2C%0A%20%20%20%20q%2C%0A%20%20%20%20rho%2C%0A%20%20%20%20sigma%2C%0A%20%20%20%20u_stack_max%2C%0A%20%20%20%20u_stack_min%2C%0A%20%20%20%20v0_samples%2C%0A%20%20%20%20x_samples%2C%0A%20%20%20%20z0_samples%2C%0A)%3A%0A%20%20%20%20from%20docs.tutorials.helper_functions.pepit_runner%20import%20PEPitRunner%0A%0A%20%20%20%20_pepit_radius_empirical%20%3D%20estimate_control_pepit_radius(%0A%20%20%20%20%20%20%20%20M%2C%0A%20%20%20%20%20%20%20%20q%2C%0A%20%20%20%20%20%20%20%20rho%2C%0A%20%20%20%20%20%20%20%20sigma%2C%0A%20%20%20%20%20%20%20%20build_lhs(P)%2C%0A%20%20%20%20%20%20%20%20u_stack_min%2C%0A%20%20%20%20%20%20%20%20u_stack_max%2C%0A%20%20%20%20%20%20%20%20z0_samples%2C%0A%20%20%20%20%20%20%20%20v0_samples%2C%0A%20%20%20%20%20%20%20%20x_samples%2C%0A%20%20%20%20)%0A%20%20%20%20_P_eigs%20%3D%20np.linalg.eigvalsh(P.toarray())%0A%20%20%20%20_mu%20%3D%20float(_P_eigs.min())%0A%20%20%20%20_L%20%3D%20float(_P_eigs.max())%0A%20%20%20%20_runner%20%3D%20PEPitRunner()%0A%20%20%20%20all_residuals_PEP_radius%2C%20_%20%3D%20_runner.osqp_admm_residuals(%0A%20%20%20%20%20%20%20%20K%3DK%2C%0A%20%20%20%20%20%20%20%20mu%3D_mu%2C%0A%20%20%20%20%20%20%20%20L%3D_L%2C%0A%20%20%20%20%20%20%20%20sigma%3Dsigma%2C%0A%20%20%20%20%20%20%20%20rho%3Drho%2C%0A%20%20%20%20%20%20%20%20R%3D_pepit_radius_empirical%2C%0A%20%20%20%20)%0A%20%20%20%20all_residuals_PEP%2C%20pepit_radius_calibrated%20%3D%20_runner.osqp_admm_residuals(%0A%20%20%20%20%20%20%20%20K%3DK%2C%0A%20%20%20%20%20%20%20%20mu%3D_mu%2C%0A%20%20%20%20%20%20%20%20L%3D_L%2C%0A%20%20%20%20%20%20%20%20sigma%3Dsigma%2C%0A%20%20%20%20%20%20%20%20rho%3Drho%2C%0A%20%20%20%20%20%20%20%20calibrate_to%3Dall_residuals%5B0%5D%2C%0A%20%20%20%20)%0A%20%20%20%20return%20all_residuals_PEP%2C%20all_residuals_PEP_radius%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22%0A%20%20%20%20%23%23%20Numerical%20Result%0A%0A%20%20%20%20We%20plot%20every%20curve%20in%20squared%20Euclidean%20units.%20The%20MIP%20curve%20uses%20the%0A%20%20%20%20dimension-dependent%20norm%20conversion%20above%3B%20simulation%20and%20PEP%20use%0A%20%20%20%20the%20squared%20Euclidean%20residual%20directly%2C%20and%20the%20simulation%20includes%20the%0A%20%20%20%20sample%20worst-case%20instance.%0A%20%20%20%20%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(%0A%20%20%20%20K%2C%0A%20%20%20%20all_residuals%2C%0A%20%20%20%20all_residuals_PEP%2C%0A%20%20%20%20all_residuals_PEP_radius%2C%0A%20%20%20%20plt%2C%0A%20%20%20%20simulation_residual_quantiles_admm%2C%0A%20%20%20%20simulation_residuals_admm%2C%0A)%3A%0A%20%20%20%20plt.figure()%0A%20%20%20%20_ks%20%3D%20range(1%2C%20K%2B1)%0A%20%20%20%20plt.plot(_ks%2C%20all_residuals_PEP%2C%20label%3D'ADMM%20(PEPit%2C%20matched%20at%20K%3D1)')%0A%20%20%20%20plt.plot(%0A%20%20%20%20%20%20%20%20_ks%2C%0A%20%20%20%20%20%20%20%20all_residuals_PEP_radius%2C%0A%20%20%20%20%20%20%20%20'%3A'%2C%0A%20%20%20%20%20%20%20%20label%3D'ADMM%20(PEPit%2C%20empirical-radius%20conditioned)'%2C%0A%20%20%20%20)%0A%20%20%20%20plt.plot(%0A%20%20%20%20%20%20%20%20_ks%2C%0A%20%20%20%20%20%20%20%20all_residuals%2C%0A%20%20%20%20%20%20%20%20label%3Dr'ADMM%20(MIP%20upper%20bound%20via%20%242n_%7Bxz%7D%5C%7C%5Ccdot%5C%7C_%5Cinfty%5E2%24)'%2C%0A%20%20%20%20)%0A%20%20%20%20_avg_line%2C%20%3D%20plt.plot(_ks%2C%20simulation_residuals_admm%2C%20'--'%2C%20label%3D'ADMM%20(in-sample%20avg)')%0A%20%20%20%20plt.fill_between(_ks%2C%20simulation_residual_quantiles_admm%5B2%5D%2C%20simulation_residual_quantiles_admm%5B0%5D%2C%20alpha%3D.18%2C%20label%3D'ADMM%20sample%2010%25-90%25%20band')%0A%20%20%20%20plt.plot(_ks%2C%20simulation_residual_quantiles_admm%5B3%5D%2C%20'%3A'%2C%20color%3D_avg_line.get_color()%2C%20label%3D'ADMM%20sample%20worst-case')%0A%20%20%20%20plt.suptitle('optimal%20control%20example')%0A%20%20%20%20plt.yscale('log')%0A%20%20%20%20plt.xlabel(r'%24K%24')%0A%20%20%20%20plt.ylabel(r'Worst-case%20squared%20fixed-point%20residual%20%24%5C%7Cr%5Ek%5C%7C_2%5E2%24')%0A%20%20%20%20plt.legend()%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