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(%0A%20%20%20%20%20%20%20%20r%22%22%22%0A%20%20%20%20We%20implement%20a%20Python%20Fusion%20model%20for%20the%20problem%20known%20as%20F-SPARC%20(fractional%20subcarrier%20and%20power%20allocation%20with%20rate%20constraints).%20Our%20implementation%20defines%20a%20mixed-integer%20exponential%20cone%20problem.%0A%0A%20%20%20%20The%20model%20comes%20from%20the%20paper%20*Bi-Perspective%20Functions%20for%20Mixed-Integer%20Fractional%20Programs%20with%20Indicator%20Variables*%20by%20Adam%20N.%20Letchford%2C%20Qiang%20N%20and%2C%20Zhaoyu%20Zhong%2C%20http%3A%2F%2Fwww.optimization-online.org%2FDB_HTML%2F2017%2F09%2F6222.html.%0A%0A%20%20%20%20Some%20test%20data%20come%20from%20http%3A%2F%2Fwww.research.lancs.ac.uk%2Fportal%2Fen%2Fdatasets%2Fofdma-optimisation(48d843b0-1fb6-4e3d-af39-ca533a37b7ab).html%0A%0A%20%20%20%20%23%20Problem%20formulation%0A%0A%20%20%20%20We%20consider%20a%20set%20%24I%24%20of%20communication%20channels%20(subcarriers)%20to%20be%20assigned%20to%20a%20set%20%24J%24%20of%20users.%20Each%20channel%20can%20be%20assigned%20to%20at%20most%20one%20user%2C%20while%20one%20user%20will%20typically%20have%20more%20channels%20assigned%20to%20them%2C%20subject%20to%20the%20requirement%20that%20user%20%24j%24%20achieves%20a%20total%20data%20rate%20of%20at%20least%20%24d_j%24%20(bits%20per%20second).%20The%20data%20rate%20of%20a%20channel%20is%20computed%20as%20follows%3A%20for%20a%20channel%20of%20bandwidth%20%24B%24%20(hertz)%2C%20noise%20level%20%24N_i%24%20(watts)%20and%20operating%20with%20power%20%24p_i%24%20(watts)%2C%20the%20maximum%20data%20rate%20is%0A%20%20%20%20%24%24f_i%20%3D%20B%20%5Clog_2(1%2B%5Cfrac%7Bp_i%7D%7BN_i%7D).%24%24%0A%0A%20%20%20%20Finally%2C%20we%20assume%20a%20minimal%20system%20power%20of%20%24%5Csigma%24%20required%20to%20keep%20the%20system%20running%20and%20a%20maximum%20power%20level%20%24P%24.%20%0A%0A%20%20%20%20The%20goal%20is%20to%20map%20the%20channels%20to%20users%20and%20assign%20their%20power%20levels%20so%20that%20the%20*energy%20efficiency*%0A%0A%20%20%20%20%24%24%5Cfrac%7B%5Ctextrm%7Btotal%20data%20rate%7D%7D%7B%5Ctextrm%7Btotal%20power%7D%7D%24%24%0A%20%20%20%20is%20maximized.%0A%0A%20%20%20%20The%20mixed-integer%20formulation%20is%20then%3A%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Brl%7D%0A%20%20%20%20%5Cmathrm%7Bmaximize%7D%20%26%20%5Cfrac%7B%5Csum_%7Bi%2Cj%7D%20B%20%5Clog_2(1%2Bp_%7Bi%2Cj%7D%2FN_i)%7D%7B%5Csigma%2B%5Csum_%7Bi%2Cj%7Dp_%7Bi%2Cj%7D%7D%20%5C%5C%0A%20%20%20%20%5Ctextrm%7Bsubject%20to%7D%20%26%20B%20%5Csum_i%20%5Clog_2(1%2Bp_%7Bi%2Cj%7D%2FN_i)%20%5Cgeq%20d_j%2C%20%5Cquad%20j%5Cin%20J%2C%20%5C%5C%0A%20%20%20%20%26%20%5Csigma%20%2B%20%5Csum_%7Bi%2Cj%7Dp_%7Bij%7D%20%5Cleq%20P%2C%20%5C%5C%0A%20%20%20%20%26%200%20%5Cleq%20p_%7Bij%7D%5Cleq%20Px_%7Bij%7D%2C%20x_%7Bij%7D%5Cin%20%5C%7B0%2C1%5C%7D%2C%20%5C%5C%0A%20%20%20%20%26%20%5Csum_j%20x_%7Bij%7D%20%5Cleq%201%2C%5Cquad%20i%5Cin%20I%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20where%20%24p_%7Bij%7D%24%20is%20the%20power%20assigned%20to%20channel%20%24i%24%20when%20used%20by%20user%20%24j%24.%20The%20first%20constraint%20describes%20the%20total%20bitrate%20of%20each%20user%2C%20the%20second%20bounds%20total%20power%2C%20and%20the%20remaining%20constraints%20with%20indicator%20variables%20model%20the%20assignment%20requirements.%0A%0A%0A%20%20%20%20%23%20Conic%20reformulation%0A%0A%20%20%20%20To%20obtain%20a%20problem%20in%20conic%20form%20we%20perform%20a%20homogenization%20as%20in%20%5BSection%205.2%20of%20Letchfort%20et%20al.%5D(https%3A%2F%2Foptimization-online.org%2F2017%2F09%2F6222%2F).%20That%20is%2C%20we%20introduce%20new%20variables%20%24t%24%2C%20%24z_%7Bij%7D%24%20and%20%24%5Ctilde%7Bp_%7Bij%7D%7D%24%20with%20the%20intention%20of%20representing%0A%0A%20%20%20%20%24%24t%20%3D%201%2F(%5Csigma%2B%5Csum_%7Bij%7Dp_%7Bij%7D)%2C%5Cquad%20z_%7Bij%7D%3DBt%5Clog_2(1%2Bp_%7Bij%7D%2FN_i)%2C%5Cquad%20%5Ctilde%7Bp_%7Bij%7D%7D%20%3D%20p_%7Bij%7Dt.%24%24%0A%0A%20%20%20%20In%20terms%20of%20the%20new%20variables%20the%20problem%20is%20equivalent%20with%0A%0A%20%20%20%20%24%24%0A%20%20%20%20%5Cbegin%7Barray%7D%7Brl%7D%0A%20%20%20%20%5Cmathrm%7Bmaximize%7D%20%26%20%5Csum_%7Bi%2Cj%7D%20z_%7Bij%7D%20%5C%5C%0A%20%20%20%20%5Ctextrm%7Bsubject%20to%7D%20%26%20z_%7Bij%7D%20%5Cleq%20B%20t%20%5Clog_2(1%2B%5Ctilde%7Bp_%7Bi%2Cj%7D%7D%2F(tN_i))%2C%20%5Cquad%20i%20%5Cin%20I%2C%20j%5Cin%20J%2C%20%5C%5C%0A%20%20%20%20%26%20d_jt%20%5Cleq%20%5Csum_i%20z_%7Bij%7D%2C%5Cquad%20j%5Cin%20J%20%5C%5C%0A%20%20%20%20%26%20t%5Csigma%20%2B%20%5Csum_%7Bi%2Cj%7D%20%5Ctilde%7Bp_%7Bi%2Cj%7D%7D%20%3D%201%2C%20%5C%5C%0A%20%20%20%20%26%201%2FP%20%5Cleq%20t%20%5Cleq%201%2F%5Csigma%2C%20%5C%5C%0A%20%20%20%20%26%200%20%5Cleq%20%5Ctilde%7Bp_%7Bi%2Cj%7D%7D%20%5Cleq%20x_%7Bij%7D%2C%20x_%7Bij%7D%5Cin%20%5C%7B0%2C1%5C%7D%2C%20%5C%5C%0A%20%20%20%20%26%20%5Csum_j%20x_%7Bij%7D%20%5Cleq%201%2C%5Cquad%20i%5Cin%20I%2C%0A%20%20%20%20%5Cend%7Barray%7D%0A%20%20%20%20%24%24%0A%0A%20%20%20%20The%20first%20constraint%20can%20equivalently%20be%20written%20as%0A%20%20%20%20%24%24t%2B%5Ctilde%7Bp_%7Bi%2Cj%7D%7D%2FN_i%20%5Cgeq%20t%20%5Cexp(%5Clog(2)z_%7Bi%2Cj%7D%2F(Bt))%24%24%0A%20%20%20%20which%2C%20using%20the%20exponential%20cone%2C%20is%0A%0A%20%20%20%20%24%24(t%2B%5Ctilde%7Bp_%7Bi%2Cj%7D%7D%2FN_i%20%2C%20t%2C%20%5Clog(2)z_%7Bi%2Cj%7D%2FB)%5Cin%20K_%5Cmathrm%7Bexp%7D.%24%24%0A%0A%20%20%20%20%23%20Implementation%20in%20Fusion%0A%0A%20%20%20%20We%20start%20that%20functions%20which%20parse%20examples%20in%20the%20format%20of%20http%3A%2F%2Fwww.research.lancs.ac.uk%2Fportal%2Fen%2Fdatasets%2Fofdma-optimisation(48d843b0-1fb6-4e3d-af39-ca533a37b7ab).html%20and%20set%20up%20some%20global%20constants.%20The%20data%20files%20contain%20user%20demands%20and%20channel%20noise%20for%20each%20example.%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%20%23%20marimo%3A%20preserve-names%0A%20%20%20%20import%20ast%2C%20sys%0A%0A%20%20%20%20constants%20%3D%20%7B%0A%20%20%20%20%20%20%22System%20Power%22%3A%2010%2C%20%0A%20%20%20%20%20%20%22Bandwidth%22%3A%201.25%2C%20%0A%20%20%20%20%20%20%22Maximum%20Power%22%3A%2036%0A%20%20%20%20%7D%0A%0A%20%20%20%20def%20parse(filename)%3A%0A%20%20%20%20%20%20%20%20data%20%3D%20%5B%5D%0A%20%20%20%20%20%20%20%20with%20open(filename%2C%20'r')%20as%20file%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20content%20%3D%20file.read()%0A%20%20%20%20%20%20%20%20%20%20%20%20if%20'Instance%3A%20'%20not%20in%20content%3A%20raise%20Exception%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20inst%20in%20content.split('Instance%3A%20')%5B1%3A%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20if%20'noise'%20not%20in%20inst%20or%20'demand'%20not%20in%20inst%3A%20raise%20Exception%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20data.append(%7B%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22noise%22%20%3A%20ast.literal_eval(inst.split('noise')%5B1%5D.split('demand')%5B0%5D.strip())%2C%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%22demand%22%20%3A%20ast.literal_eval(inst.split('demand')%5B1%5D.strip())%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%7D)%0A%20%20%20%20%20%20%20%20return%20data%20%0A%0A%20%20%20%20BW%20%3D%20constants%5B%22Bandwidth%22%5D%0A%20%20%20%20P%20%3D%20constants%5B%22Maximum%20Power%22%5D%0A%20%20%20%20SIGMA%20%3D%20constants%5B%22System%20Power%22%5D%0A%20%20%20%20return%20BW%2C%20P%2C%20SIGMA%2C%20parse%2C%20sys%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22The%20next%20function%20is%20a%20direct%20implementation%20of%20the%20conic%20model%20above.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(BW%2C%20P%2C%20SIGMA)%3A%0A%20%20%20%20from%20mosek.fusion%20import%20Model%2C%20Domain%2C%20Expr%2C%20ObjectiveSense%2C%20Matrix%2C%20Var%2C%20SolutionStatus%2C%20%20AccSolutionStatus%0A%20%20%20%20import%20mosek.fusion.pythonic%0A%20%20%20%20from%20math%20import%20log%0A%0A%20%20%20%20def%20fsparc_model(n%2C%20d)%3A%0A%20%20%20%20%20%20%20%20M%20%3D%20Model()%0A%20%20%20%20%20%20%20%20I%2C%20J%20%3D%20(len(n)%2C%20len(d))%0A%20%20%20%20%20%20%20%20z%20%3D%20M.variable('z'%2C%20%5BI%2C%20J%5D%2C%20Domain.greaterThan(0.0))%0A%20%20%20%20%20%20%20%20x%20%3D%20M.variable('x'%2C%20%5BI%2C%20J%5D%2C%20Domain.binary())%0A%20%20%20%20%20%20%20%20p_tilde%20%3D%20M.variable('p_tilde'%2C%20%5BI%2C%20J%5D%2C%20Domain.greaterThan(0.0))%0A%20%20%20%20%20%20%20%20t%20%3D%20M.variable('t'%2C%201%2C%20Domain.inRange(1.0%20%2F%20P%2C%201.0%20%2F%20SIGMA))%0A%20%20%20%20%20%20%20%20M.constraint(Expr.sum(x%2C%201)%20%3C%3D%201.0)%0A%20%20%20%20%20%20%20%20M.constraint(SIGMA%20*%20t%20%2B%20Expr.sum(p_tilde)%20%3D%3D%201.0)%0A%20%20%20%20%20%20%20%20M.constraint(p_tilde%20%3C%3D%20x)%0A%20%20%20%20%20%20%20%20M.constraint(Expr.sum(z%2C%200)%20%3E%3D%20Expr.mulElm(Expr.repeat(t%2C%20J%2C%200)%2C%20d))%0A%20%20%20%20%20%20%20%20M.objective('obj'%2C%20ObjectiveSense.Maximize%2C%20Expr.sum(z))%0A%20%20%20%20%20%20%20%20p_over_n%20%3D%20Expr.mulElm(p_tilde%2C%20%5B%5B1%20%2F%20n_i%5D%20*%20J%20for%20n_i%20in%20n%5D)%0A%20%20%20%20%20%20%20%20multi_t%20%3D%20Expr.repeat(Expr.repeat(t%2C%20I%2C%200)%2C%20J%2C%201)%0A%20%20%20%20%20%20%20%20M.constraint(Expr.stack(2%2C%20multi_t%20%2B%20p_over_n%2C%20multi_t%2C%20z%20*%20(log(2.0)%20%2F%20BW))%2C%20Domain.inPExpCone())%0A%20%20%20%20%20%20%20%20return%20M%0A%20%20%20%20return%20AccSolutionStatus%2C%20fsparc_model%2C%20log%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22Next%20we%20run%20the%20model%20and%20present%20some%20statistics%20for%20a%20few%20examples.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(BW%2C%20P%2C%20SIGMA%2C%20log)%3A%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%0A%20%20%20%20def%20stats(M%2C%20n%2C%20d)%3A%0A%20%20%20%20%20%20%20%20I%2C%20J%20%3D%20(len(n)%2C%20len(d))%0A%20%20%20%20%20%20%20%20t%20%3D%20M.getVariable('t').level()%5B0%5D%0A%20%20%20%20%20%20%20%20pval%20%3D%20M.getVariable('p_tilde').level().reshape(%5BI%2C%20J%5D)%20%2F%20t%0A%20%20%20%20%20%20%20%20xval%20%3D%20M.getVariable('x').level().reshape(%5BI%2C%20J%5D)%0A%20%20%20%20%20%20%20%20power%20%3D%20%5Bsum(pval%5Bi%2C%20%3A%5D)%20for%20i%20in%20range(I)%5D%0A%20%20%20%20%20%20%20%20user%20%3D%20%5Bint(sum(%5Bxval%5Bi%5D%5Bj%5D%20*%20j%20for%20j%20in%20range(J)%5D))%20for%20i%20in%20range(I)%5D%0A%20%20%20%20%20%20%20%20print('Total%20data%20rate%3A%20demanded%20%7B0%7D%2C%20provided%20%7B1%7D'.format(sum(d)%2C%20sum(%5BBW%20*%20log(1%20%2B%20power%5Bi%5D%20%2F%20n%5Bi%5D)%20%2F%20log(2.0)%20for%20i%20in%20range(I)%5D)))%0A%20%20%20%20%20%20%20%20print('Power%3A%20system%20%7B0%3A.3f%7D%2C%20used%20%7B1%3A.3f%7D%2C%20maximal%3A%20%7B2%3A.3f%7D'.format(SIGMA%2C%20SIGMA%20%2B%20sum(power)%2C%20P))%0A%20%20%20%20%20%20%20%20colors%20%3D%20%5B'red'%2C%20'green'%2C%20'blue'%2C%20'cyan'%2C%20'yellow'%2C%20'violet'%5D%0A%20%20%20%20%20%20%20%20plt.bar(range(I)%2C%20power%2C%20color%3D%5Bcolors%5Buser%5Bi%5D%5D%20for%20i%20in%20range(I)%5D)%0A%20%20%20%20return%20(stats%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%20Example%201%0A%0A%20%20%20%20We%20run%20an%20example%20with%20realistic%20data%20(72%20channels%2C%204%20users)%20with%20a%20very%20short%20time%20limit.%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_(AccSolutionStatus%2C%20fsparc_model%2C%20parse%2C%20stats%2C%20sys)%3A%0A%20%20%20%20def%20_()%3A%0A%20%20%20%20%20%20%20%20_instance%20%3D%20parse('data%2Flancaster%2F4_0.8.txt')%5B0%5D%0A%20%20%20%20%20%20%20%20n%2C%20d%20%3D%20(_instance%5B'noise'%5D%2C%20_instance%5B'demand'%5D)%0A%20%20%20%20%20%20%20%20M%20%3D%20fsparc_model(n%2C%20d)%0A%20%20%20%20%20%20%20%20M.setLogHandler(sys.stdout)%0A%20%20%20%20%20%20%20%20M.setSolverParam('optimizerMaxTime'%2C%2040.0)%0A%20%20%20%20%20%20%20%20M.solve()%0A%20%20%20%20%20%20%20%20M.acceptedSolutionStatus(AccSolutionStatus.Feasible)%0A%0A%20%20%20%20%20%20%20%20print(M.getPrimalSolutionStatus())%0A%20%20%20%20%20%20%20%20return%20stats(M%2C%20n%2C%20d)%0A%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%20The%20bar%20plot%20shows%20power%20assignments%20to%20the%20channels.%20Each%20color%20represents%20one%20user.%0A%0A%20%20%20%20%23%20Example%202%0A%0A%20%20%20%20The%20next%20example%20is%20a%20small%2C%20randomly%20generated%20set%20of%20data%20with%2015%20channels%20and%203%20users%2C%20which%20we%20should%20solve%20to%20optimality.%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_(fsparc_model%2C%20parse%2C%20stats%2C%20sys)%3A%0A%20%20%20%20_instance%20%3D%20parse('data%2Fsmall-random%2Frandom_15_3_0.85.txt')%5B0%5D%0A%20%20%20%20n%2C%20d%20%3D%20(_instance%5B'noise'%5D%2C%20_instance%5B'demand'%5D)%0A%20%20%20%20M%20%3D%20fsparc_model(n%2C%20d)%0A%20%20%20%20M.setLogHandler(sys.stdout)%0A%20%20%20%20M.solve()%0A%20%20%20%20stats(M%2C%20n%2C%20d)%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%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%0Aif%20__name__%20%3D%3D%20%22__main__%22%3A%0A%20%20%20%20app.run()%0A
3eea51312f2ff31370dc1207436267da03f9b426c2ca0dfcde0350ec81eebdee