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%20%23%23%20The%20Unit%20Commitment%20Problem%20(UCP)%0A%0A%20%20%20%20The%20Unit%20Commitment%20Problem%20is%20an%20optimization%20problem%20in%20electrical%20power%20production.%20The%20goal%20is%20to%20assign%20production%20levels%20and%20switch-on%2Fswitch-off%20times%20to%20generators%20in%20a%20network%20so%20that%20the%20total%20production%20at%20any%20time%20meets%20the%20predicted%20demand%20for%20electricity.%20In%20this%20demonstration%20we%20consider%20minimizing%20production%20cost%20subject%20to%20constraints%20on%20the%20generating%20units%3A%20minimal%2Fmaximal%20production%20levels%2C%20minimum%20uptime%2Fdowntime%20and%20maximum%20ramp-up%2Framp-down.%20The%20cost%20consists%20of%20quadratic%20production%20costs%2C%20fixed%20operating%20cost%20and%20start-up%20costs.%20That%20makes%20our%20version%20of%20UCP%20a%20Mixed%20Integer%20Conic%20Quadratic%20Optimization%20problem%20(MICQO).%0A%0A%20%20%20%20These%20are%20the%20most%20basic%20constraint%20types%20for%20a%20UCP%20model%2C%20see%20for%20example%20%5Bthis%20paper%5D(http%3A%2F%2Fpierrepinson.com%2Fdocs%2FHeideJorgensenetal2016.pdf)%20and%20references%20therein.%20We%20only%20model%20the%20total%20power%20generation%20of%20thermal%20units%20with%20constant%20startup%20costs%2C%20but%20not%20the%20storage%20capacity%20of%20hydro%20units%20nor%20cascades%2C%20power%20flow%20or%20bidding%20aspects%3B%20see%20for%20instance%20%5Ba%20more%20thorough%20real-world%20model%5D(https%3A%2F%2Fwww.ferc.gov%2Flegal%2Fstaff-reports%2Frto-COMMITMENT-TEST.pdf).%20Many%20of%20those%20features%20could%20be%20added%20by%20extending%20our%20linear%20model.%20We%20use%20examples%20derived%20from%20datasets%20from%20the%20%5BOR-library%5D(http%3A%2F%2Fpeople.brunel.ac.uk%2F~mastjjb%2Fjeb%2Finfo.html)%20(data%20is%20modified%20to%20fit%20into%20our%20framework).%20They%20are%20loaded%20with%20the%20script%3A%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%20ucpParser%0A%20%20%20%20return%20(ucpParser%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%20The%20source%20file%20and%20examples%20are%20available%20in%20the%20%5BGitHub%5D(https%3A%2F%2Fgithub.com%2FMOSEK%2FTutorials%2Ftree%2Fmaster%2Funitcommitment)%20repository%20containing%20this%20notebook.%20%0A%0A%20%20%20%20We%20demonstrate%20such%20aspects%20of%20working%20with%20MOSEK%20as%3A%0A%0A%20%20%20%20*%20setting%20up%20a%20model%2C%0A%20%20%20%20*%20using%20variable%20slices%20to%20effectively%20write%20constraints%2C%0A%20%20%20%20*%20retrieving%20optimizer%20statistics%2C%0A%20%20%20%20*%20reoptimization.%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%20'%25matplotlib%20inline'%20command%20supported%20automatically%20in%20marimo%0A%20%20%20%20import%20matplotlib.pyplot%20as%20plt%0A%20%20%20%20import%20matplotlib.cm%20as%20cm%0A%20%20%20%20import%20numpy%20as%20np%0A%20%20%20%20import%20sys%0A%20%20%20%20from%20mosek.fusion%20import%20Model%2C%20Domain%2C%20Expr%2C%20ObjectiveSense%2C%20Matrix%2C%20Var%2C%20SolutionType%2C%20ProblemStatus%2C%20AccSolutionStatus%0A%20%20%20%20import%20mosek.fusion.pythonic%0A%20%20%20%20return%20(%0A%20%20%20%20%20%20%20%20AccSolutionStatus%2C%0A%20%20%20%20%20%20%20%20Domain%2C%0A%20%20%20%20%20%20%20%20Expr%2C%0A%20%20%20%20%20%20%20%20Model%2C%0A%20%20%20%20%20%20%20%20ObjectiveSense%2C%0A%20%20%20%20%20%20%20%20ProblemStatus%2C%0A%20%20%20%20%20%20%20%20SolutionType%2C%0A%20%20%20%20%20%20%20%20Var%2C%0A%20%20%20%20%20%20%20%20cm%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%20sys%2C%0A%20%20%20%20)%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%23%20Problem%20formulation%0A%0A%20%20%20%20In%20the%20UCP%20problem%20we%20are%20given%20a%20number%20%24N%24%20of%20generators%20(*units*)%2C%20a%20time%20horizon%20of%20%24T%24%20*periods*%20(usually%20hours%20or%20half-hours)%20and%20initial%20conditions%20(see%20later).%20The%20goal%20is%20to%20determine%20electricity%20production%20levels%20for%20all%20units%20and%20all%20periods%20so%20that%20a%20total%20predicted%20*demand*%20is%20satisfied%20in%20each%20period.%20The%20operation%20of%20the%20units%20is%20subject%20to%20technical%20constraints%20which%20will%20be%20introduced%20below.%20The%20operating%20costs%20are%20known%20for%20each%20unit%20and%20can%20be%20split%20into%20*production%20costs*%2C%20*fixed%20operating%20costs*%20when%20the%20unit%20is%20on%20and%20*startup%20costs*%20incurred%20when%20the%20unit%20is%20switched%20on.%0A%0A%20%20%20%20We%20can%20model%20this%20setup%20using%20these%20simple%20variables%3A%0A%0A%20%20%20%20*%20a%20continous%20variable%20%24p%5Cin%5Cmathbb%7BR%7D_%2B%5E%7BN%5Ctimes%20T%7D%24%20such%20that%20%24p_%7Bi%2Ct%7D%24%20is%20the%20power%20generation%20of%20unit%20%24i%24%20at%20time%20%24t%24%2C%0A%20%20%20%20*%20a%20binary%20variable%20%24u%5Cin%5C%7B0%2C1%5C%7D%5E%7BN%5Ctimes%20T%7D%24%2C%20where%20%24u_%7Bi%2Ct%7D%24%20indicates%20if%20unit%20%24i%24%20is%20*on*%20(producing%20energy)%20at%20time%20%24t%24%2C%0A%20%20%20%20*%20a%20binary%20variable%20%24v%5Cin%5C%7B0%2C1%5C%7D%5E%7BN%5Ctimes%20T%7D%24%2C%20where%20%24v_%7Bi%2Ct%7D%24%20indicates%20if%20unit%20%24i%24%20has%20been%20started%20(switched%20from%20*off*%20to%20*on*)%20at%20time%20%24t%24%2C%20in%20other%20words%20%24v_%7Bi%2Ct%7D%3D%5Cmax(u_%7Bi%2Ct%7D-u_%7Bi-1%2Ct%7D%2C%5C%200)%24.%0A%0A%20%20%20%20It%20will%20be%20convenient%20to%20declare%20variables%20of%20size%20%24N%5Ctimes(T%2B1)%24%20and%20use%20the%20first%20column%20for%20initial%20conditions%2C%20that%20is%20production%20levels%20and%20on%2Foff%20status%20at%20the%20beginning%20of%20the%20simulation.%20The%20code%20that%20initializes%20a%20Fusion%20model%20is%20shown%20below.%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%20Model)%3A%0A%20%20%20%20class%20UCP%3A%0A%20%20%20%20%20%20%20%20%23%20Initialize%20with%20input%20data%0A%20%20%20%20%20%20%20%20def%20__init__(self%2C%20T%2C%20N)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20self.T%2C%20self.N%20%3D%20T%2C%20N%0A%20%20%20%20%20%20%20%20%20%20%20%20self.M%20%3D%20Model()%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Production%20level%20---%20p_ih%0A%20%20%20%20%20%20%20%20%20%20%20%20self.p%20%3D%20self.M.variable(%5Bself.N%2C%20self.T%2B1%5D%2C%20Domain.greaterThan(0.0))%20%20%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Is%20active%3F%20---%20u_ih%0A%20%20%20%20%20%20%20%20%20%20%20%20self.u%20%3D%20self.M.variable(%5Bself.N%2C%20self.T%2B1%5D%2C%20Domain.binary())%20%20%20%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Is%20at%20startup%3F%20---%20v_ih%20%3E%3D%20u_ih-u_(i-1)h%0A%20%20%20%20%20%20%20%20%20%20%20%20self.v%20%3D%20self.M.variable(%5Bself.N%2C%20self.T%2B1%5D%2C%20Domain.binary())%20%20%20%0A%20%20%20%20%20%20%20%20%20%20%20%20self.M.constraint(self.v%5B%3A%2C1%3A%5D%20%3E%3D%20self.u%5B%3A%2C1%3A%5D%20-%20self.u%5B%3A%2C%3A-1%5D)%0A%20%20%20%20return%20(UCP%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22Note%20how%20we%20use%20slicing%20to%20express%20the%20constraint%20%24v_%7Bi%2Ct%7D%5Cgeq%20u_%7Bi%2Ct%7D-u_%7Bi-1%2Ct%7D%24%20simultaneously%20for%20all%20%24i%2Ct%24.%20The%20two%20slices%20of%20%24u%24%20have%20the%20same%20dimensions%20but%20are%20shifted%20by%20%241%24%20with%20respect%20to%20each%20other.%20The%20%24i%2Ct%24-entry%20of%20the%20constraint%20expression%20is%20now%20just%20%24v_%7Bi%2Ct%7D-(u_%7Bi%2Ct%7D-u_%7Bi-1%2Ct%7D)%5Cgeq%200%24%20as%20required.%20We%20are%20going%20to%20use%20this%20type%20of%20slicing%20heavily%20in%20what%20follows.%20Beacuse%20most%20of%20the%20operational%20constraints%20of%20the%20units%20will%20need%20to%20be%20broadcasted%20to%20all%20time%20periods%2C%20this%20repeat%20method%20will%20come%20in%20handy.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(np)%3A%0A%20%20%20%20%23%20Create%20a%20matrix%20with%20R%20columns%20each%20equal%20to%20v%20(horizontal%20repeat)%0A%20%20%20%20def%20repeat(v%2C%20R)%3A%0A%20%20%20%20%20%20%20%20return%20np.repeat(v%2C%20R).reshape(len(v)%2C%20R)%0A%20%20%20%20return%20(repeat%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%23%23%20Production%20constraints%0A%0A%20%20%20%20Each%20unit%2C%20when%20in%20use%2C%20operates%20between%20a%20minimal%20and%20maximal%20power%20level%20%24p%5E%7B%5Ctextrm%7Bmin%7D%7D_i%24%20and%20%24p%5E%7B%5Ctextrm%7Bmax%7D%7D_i%24.%20If%20the%20unit%20is%20off%2C%20its%20production%20is%20%240%24%3A%0A%0A%20%20%20%20%24%24p%5E%7B%5Ctextrm%7Bmin%7D%7D_i%20u_%7Bi%2Ct%7D%5Cleq%20p_%7Bi%2Ct%7D%5Cleq%20p%5E%7B%5Ctextrm%7Bmax%7D%7D_i%20u_%7Bi%2Ct%7D.%24%24%0A%0A%20%20%20%20Moreover%2C%20at%20each%20time%20%24t%24%20the%20total%20production%20must%20satisfy%20the%20current%20demand%20%24d_t%24%3A%0A%0A%20%20%20%20%24%24%5Csum_i%20p_%7Bi%2Ct%7D%20%5Cgeq%20d_t%2C%5Cquad%20t%3D1%2C%5Cldots%2CT.%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%0Adef%20_(Domain%2C%20Expr%2C%20repeat)%3A%0A%20%20%20%20%23%20Production%20constraints%0A%20%20%20%20def%20consProd(ucp%2C%20pmin%2C%20pmax%2C%20demand)%3A%0A%20%20%20%20%20%20%20%20N%2C%20T%2C%20p%2C%20u%2C%20v%20%3D%20ucp.N%2C%20ucp.T%2C%20ucp.p%2C%20ucp.u%2C%20ucp.v%0A%20%20%20%20%20%20%20%20%23%20Maximal%20and%20minimal%20production%20of%20each%20unit%0A%20%20%20%20%20%20%20%20ucp.M.constraint(p%20%3C%3D%20Expr.mulElm(repeat(pmax%2C%20T%2B1)%2C%20u))%0A%20%20%20%20%20%20%20%20ucp.M.constraint(p%20%3E%3D%20Expr.mulElm(repeat(pmin%2C%20T%2B1)%2C%20u))%0A%20%20%20%20%20%20%20%20%23%20Total%20demand%20in%20periods%20is%20achieved%0A%20%20%20%20%20%20%20%20ucp.M.constraint(Expr.sum(p%2C%200)%5B1%3A%5D%2C%20Domain.greaterThan(demand))%0A%20%20%20%20return%20(consProd%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22Here%20we%20used%20the%20%60%60Expr.sum(p%2C%200)%60%60%20operator%20to%20sum%20the%20%24p%24%20matrix%20along%20the%20%240%24-th%20dimension%20i.e.%20to%20compute%20sums%20within%20each%20column.%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%20%23%23%23%20Ramp%20constraints%0A%0A%20%20%20%20Ramp-up%20and%20ramp-down%20constraints%20indicate%20that%20the%20output%20of%20a%20generator%20cannot%20vary%20too%20rapidly%20between%20periods.%20%0A%0A%20%20%20%20%24%24p_%7Bi-1%2Ct%7D-r%5E%7B%5Ctextrm%7Bdown%7D%7D_i%5Cleq%20p_%7Bi%2Ct%7D%5Cleq%20p_%7Bi-1%2Ct%7D%2Br%5E%7B%5Ctextrm%7Bup%7D%7D_i.%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%0Adef%20_(repeat)%3A%0A%20%20%20%20%23%20Ramp%20constraints%20%0A%20%20%20%20def%20consRamp(ucp%2C%20rdown%2C%20rup)%3A%0A%20%20%20%20%20%20%20%20N%2C%20T%2C%20p%2C%20u%2C%20v%20%3D%20ucp.N%2C%20ucp.T%2C%20ucp.p%2C%20ucp.u%2C%20ucp.v%0A%20%20%20%20%20%20%20%20Delta%20%3D%20p%5B%3A%2C1%3A%5D%20-%20p%5B%3A%2C%3A-1%5D%0A%20%20%20%20%20%20%20%20ucp.M.constraint(Delta%20%3C%3D%20repeat(rup%2C%20T))%0A%20%20%20%20%20%20%20%20ucp.M.constraint(Delta%20%3E%3D%20-repeat(rdown%2C%20T))%0A%20%20%20%20return%20(consRamp%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%23%23%20Initial%20values%0A%0A%20%20%20%20We%20also%20input%20the%20state%20of%20the%20system%20in%20the%20last%20period%20before%20the%20start%20of%20the%20simulation.%20We%20will%20load%20that%20information%20by%20fixing%20the%20values%20in%20the%20%240%24-th%20column%20of%20the%20variables.%20For%20each%20generator%20the%20initial%20state%20is%20described%20by%20the%20values%20%24p_%7Bi%2C0%7D%24%20and%20%24u_%7Bi%2C0%7D%24%20as%20well%20as%20%24l_i%24%20-%20the%20number%20of%20periods%20the%20unit%20has%20already%20spent%20in%20its%20current%20state%20(on%2Foff).%20This%20will%20be%20required%20to%20properly%20handle%20the%20initial%20uptime%2Fdowntime%20constraint%20(see%20next).%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_(Var)%3A%0A%20%20%20%20%23%20Initial%20values%0A%20%20%20%20def%20consInit(ucp%2C%20p0%2C%20u0%2C%20l0)%3A%0A%20%20%20%20%20%20%20%20N%2C%20T%2C%20p%2C%20u%2C%20v%20%3D%20ucp.N%2C%20ucp.T%2C%20ucp.p%2C%20ucp.u%2C%20ucp.v%0A%20%20%20%20%20%20%20%20%23%20Fix%20production%20in%20the%20immediately%20preceeding%20period%0A%20%20%20%20%20%20%20%20ucp.M.constraint(Var.flatten(p%5B%3A%2C0%5D)%20%3D%3D%20p0)%0A%20%20%20%20%20%20%20%20ucp.M.constraint(Var.flatten(u%5B%3A%2C0%5D)%20%3D%3D%20u0)%0A%20%20%20%20return%20(consInit%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%23%23%20Minimal%20downtime%2Fuptime%0A%0A%20%20%20%20When%20a%20large%20thermal%20unit%20is%20turned%20on%2Foff%20it%20will%20often%20be%20impossible%20to%20turn%20it%20back%20off%2Fon%20before%20a%20minimal%20uptime%2Fdowntime%20%24m%5E%7B%5Ctextrm%7Bup%7D%7D_i%2Fm%5E%7B%5Ctextrm%7Bdown%7D%7D_i%24%20(measured%20in%20periods)%20has%20elapsed.%20%0A%0A%20%20%20%20These%20conditions%20are%20modelled%20with%20the%20binary%20variable%20%24u%24.%20If%20the%20%24i%24-th%20unit%20is%20off%20at%20two%20times%20%24t%24%20and%20%24t%2Bw%24%20with%20%24w%3Cm%5E%7B%5Ctextrm%7Bup%7D%7D_i%24%20then%20it%20must%20be%20off%20at%20all%20times%20%24s%24%20in%20between%20not%20to%20violate%20the%20minimum%20uptime%20condition.%20This%20is%20equivalent%20to%0A%0A%20%20%20%20%24%24u_%7Bi%2Ct%7D%2Bu_%7Bi%2Ct%2Bw%7D%5Cgeq%20u_%7Bi%2Ct%2Bs%7D%2C%5Cquad%201%20%3C%20s%3C%20w%20%3C%20m%5E%7B%5Ctextrm%7Bup%7D%7D_i.%24%24%0A%0A%20%20%20%20This%20condition%20can%20again%20be%20expressed%20by%20combining%20three%20slices%20of%20%24u%24%20with%20offsets%20%240%2Cs%2Cw%24%20(this%20time%20separately%20for%20each%20row%20of%20%24u%24).%20A%20similar%20linear%20inequality%20is%20used%20to%20model%20minimum%20downtime.%0A%0A%20%20%20%20Finally%20the%20initial%20condition%20%24l_i%24%20may%20require%20fixing%20some%20%24u_%7Bi%2Ct%7D%24%20for%20small%20%24t%24%20if%20historically%20a%20unit%20was%20not%20on%2Foff%20long%20enough%20to%20fill%20the%20minimum%20uptime%2Fdowntime%20conditions%20at%20the%20beginning%20of%20the%20process.%0A%20%20%20%20%22%22%22%0A%20%20%20%20)%0A%20%20%20%20return%0A%0A%0A%40app.function%0A%23%20Minimal%20downtime%2Fuptime%20requirements%0Adef%20consMinT(ucp%2C%20mdown%2C%20mup%2C%20u0%2C%20l0)%3A%0A%20%20%20%20N%2C%20T%2C%20p%2C%20u%2C%20v%20%3D%20ucp.N%2C%20ucp.T%2C%20ucp.p%2C%20ucp.u%2C%20ucp.v%0A%20%20%20%20%23%20Different%20units%20will%20have%20different%20requirements%0A%20%20%20%20for%20unit%20in%20range(N)%3A%0A%20%20%20%20%20%20%20%20%23%20First%2C%20take%20care%20of%20neighbourhood%20of%20initial%20conditions%0A%20%20%20%20%20%20%20%20if%20u0%5Bunit%5D%20%3D%3D%200%20and%20l0%5Bunit%5D%20%3C%20mdown%5Bunit%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20ucp.M.constraint(u%5Bunit%2C1%3Amdown%5Bunit%5D-l0%5Bunit%5D%2B1%5D%20%3D%3D%200)%0A%20%20%20%20%20%20%20%20if%20u0%5Bunit%5D%20%3D%3D%201%20and%20l0%5Bunit%5D%20%3C%20mup%5Bunit%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20ucp.M.constraint(u%5Bunit%2C1%3Amup%5Bunit%5D-l0%5Bunit%5D%2B1%5D%20%3D%3D%201)%0A%20%20%20%20%20%20%20%20%23%20If%20we%20are%20off%20in%20i%20and%20i%2Bw%20(w%3Cmup)%2C%20we%20must%20be%20off%20in%20between%3A%0A%20%20%20%20%20%20%20%20for%20w%20in%20range(1%2Cmup%5Bunit%5D)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20s%20in%20range(1%2Cw)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ucp.M.constraint(u%5Bunit%2C0%3AT-w%2B1%5D%20%2B%20u%5Bunit%2Cw%3AT%2B1%5D%20%3E%3D%20u%5Bunit%2Cs%3AT-w%2B1%2Bs%5D)%0A%20%20%20%20%20%20%20%20%23%20If%20we%20are%20on%20in%20i%20and%20i%2Bw%20(w%3Cmdown)%2C%20we%20must%20be%20on%20in%20between%3A%0A%20%20%20%20%20%20%20%20for%20w%20in%20range(1%2Cmdown%5Bunit%5D)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20for%20s%20in%20range(1%2Cw)%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20ucp.M.constraint(u%5Bunit%2C0%3AT-w%2B1%5D%20%2B%20u%5Bunit%2Cw%3AT%2B1%5D%20%3C%3D%20u%5Bunit%2Cs%3AT-w%2B1%2Bs%5D%20%2B%201)%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%23%23%20Objective%3A%20cost%20minimization%0A%0A%20%20%20%20We%20consider%20three%20types%20of%20power%20generation%20costs%20per%20unit%20%24i%24%20per%20period%3A%0A%0A%20%20%20%20*%20production%20cost%20%24C_i(p)%20%3D%20a_ip%5E2%2Bb_ip%2Bc_i%24%20depending%20quadratically%20on%20power%20generation%20%24p%24%2C%0A%20%20%20%20*%20a%20fixed%20cost%20of%20%24%5Ctextrm%7Bfx%7D_i%24%20when%20in%20production%20(%24u_%7Bi%2Ct%7D%3D1%24)%2C%0A%20%20%20%20*%20a%20fixed%20startup%20cost%20%24%5Ctextrm%7Bsc%7D_i%24%20whenever%20the%20unit%20is%20being%20switched%20on%20(%24v_%7Bi%2Ct%7D%3D1%24).%0A%0A%20%20%20%20The%20total%20cost%20is%20the%20sum%20of%20all%20these%20costs.%20The%20costs%20associated%20to%20%24b%2C%20%5Ctextrm%7Bfx%7D%2C%20%5Ctextrm%7Bsc%7D%24%20can%20be%20added%20up%20simply%20using%20inner%20product%20(note%20that%20the%20%60%60Expr.dot%60%60%20does%20the%20right%20thing%20for%20any%20two%20objects%20of%20the%20same%20shape).%20The%20quadratic%20terms%20(optional)%20are%20upper-bounded%20using%20auxiliary%20variables%20%24t_i%24%20with%20rotated%20quadratic%20cone%20inequalities%3A%0A%0A%20%20%20%20%24%242%5Ccdot%5Cfrac12%5Ccdot%20t_i%5Cgeq%20%5Csum_%7Bt%7Dp_%7Bi%2Ct%7D%5E2.%24%24%0A%0A%20%20%20%20There%20is%20one%20such%20constraint%20per%20row%2C%20but%20they%20are%20bundled%20together%20in%20matrix%20form.%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%20ObjectiveSense%2C%20np%2C%20repeat)%3A%0A%20%20%20%20%23%20Objective%0A%20%20%20%20def%20objective(ucp%2C%20a%2C%20b%2C%20c%2C%20fx%2C%20sc%2C%20onlyLinear%3DFalse)%3A%0A%20%20%20%20%20%20%20%20N%2C%20T%2C%20p%2C%20u%2C%20v%20%3D%20ucp.N%2C%20ucp.T%2C%20ucp.p%2C%20ucp.u%2C%20ucp.v%0A%20%20%20%20%20%20%20%20%23%20Startup%20costs%20and%20fixed%20operating%20costs%0A%20%20%20%20%20%20%20%20scCost%20%3D%20Expr.dot(v%5B%3A%2C1%3A%5D%2C%20repeat(sc%2C%20T))%0A%20%20%20%20%20%20%20%20fxCost%20%3D%20Expr.dot(u%5B%3A%2C1%3A%5D%2C%20repeat(fx%2C%20T))%0A%20%20%20%20%20%20%20%20%23%20Linear%20part%20of%20production%20cost%20bp%0A%20%20%20%20%20%20%20%20b%20%3D%20np.reshape(np.array(b)%2C%20(len(b)%2C%201))%0A%20%20%20%20%20%20%20%20bCost%20%3D%20Expr.sum(b.T%20%40%20p%5B%3A%2C1%3A%5D)%0A%20%20%20%20%20%20%20%20allCosts%20%3D%20%5BbCost%2C%20Expr.constTerm(sum(c)*ucp.T)%2C%20fxCost%2C%20scCost%5D%0A%20%20%20%20%20%20%20%20%23%20Quadratic%20part%20of%20production%20cost%20ap%5E2%0A%20%20%20%20%20%20%20%20if%20not%20onlyLinear%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20t%20%3D%20ucp.M.variable(ucp.N)%20%0A%20%20%20%20%20%20%20%20%20%20%20%20ucp.M.constraint(Expr.hstack(Expr.constTerm(N%2C%200.5)%2C%20t%2C%20p%5B%3A%2C1%3A%5D)%2C%20Domain.inRotatedQCone())%0A%20%20%20%20%20%20%20%20%20%20%20%20aCost%20%3D%20Expr.dot(a%2C%20t)%0A%20%20%20%20%20%20%20%20%20%20%20%20allCosts.append(aCost)%0A%20%20%20%20%20%20%20%20%23%20Total%0A%20%20%20%20%20%20%20%20ucp.M.objective(ObjectiveSense.Minimize%2C%20Expr.add(allCosts))%0A%20%20%20%20return%20(objective%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%23%20Example%0A%0A%20%20%20%20The%20available%20datasets%20span%20over%20a%20time%20horizon%20of%20%24T%3D24%24%20hours%2C%20but%20the%20parser%20allows%20to%20extend%20them%20by%20repeating%20periodically%20over%20a%20given%20number%20of%20days.%20We%20load%20the%20data%20and%20apply%20all%20the%20constraints%20defined%20previously.%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_(UCP%2C%20consInit%2C%20consProd%2C%20consRamp%2C%20objective%2C%20ucpParser)%3A%0A%20%20%20%20example%3D'RCUC%2FT-Ramp%2F20_0_3_w.mod'%0A%20%20%20%20numDays%20%3D%203%0A%0A%20%20%20%20T%2C%20N%2C%20demand%2C%20pmin%2C%20pmax%2C%20rdown%2C%20rup%2C%20mdown%2C%20mup%2C%20a%2C%20b%2C%20c%2C%20sc%2C%20fx%2C%20p0%2C%20u0%2C%20l0%20%3D%20ucpParser.loadModel(example%2C%20numDays)%0A%0A%20%20%20%20ucp%20%3D%20UCP(T%2C%20N)%0A%20%20%20%20consProd(ucp%2C%20pmin%2C%20pmax%2C%20demand)%0A%20%20%20%20consRamp(ucp%2C%20rdown%2C%20rup)%0A%20%20%20%20consInit(ucp%2C%20p0%2C%20u0%2C%20l0)%0A%20%20%20%20consMinT(ucp%2C%20mdown%2C%20mup%2C%20u0%2C%20l0)%0A%20%20%20%20objective(ucp%2C%20a%2C%20b%2C%20c%2C%20fx%2C%20sc)%0A%20%20%20%20return%20demand%2C%20pmax%2C%20ucp%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22We%20can%20now%20solve%20the%20model.%20It%20is%20recommended%20to%20relax%20termination%20critieria%20for%20the%20mixed-integer%20optimizer%20using%20MOSEK%20parameters.%20Here%20we%20will%20accept%20solutions%20with%20relative%20optimality%20tolerance%20%240.6%5C%25%24%2C%20i.e.%20%240.6%5C%25%24%20from%20the%20objective%20bound.%20Another%20option%20is%20to%20set%20a%20time%20limit.%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(sys%2C%20ucp)%3A%0A%20%20%20%20ucp.M.setLogHandler(sys.stdout)%0A%20%20%20%20ucp.M.setSolverParam(%22mioTolRelGap%22%2C%200.03)%20%23%203%25%0A%20%20%20%20%23ucp.M.setSolverParam(%22mioMaxTime%22%2C%2030.0)%20%20%20%23%2030%20sec.%0A%20%20%20%20ucp.M.solve()%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%22Next%2C%20we%20print%20some%20solution%20statistics%20and%20fetch%20the%20solution%20values%20for%20%24p%24%20and%20%24u%24%3A%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(AccSolutionStatus%2C%20ProblemStatus%2C%20SolutionType%2C%20ucp)%3A%0A%20%20%20%20%23%20Fetch%20results%0A%20%20%20%20def%20results(ucp)%3A%0A%20%20%20%20%20%20%20%20N%2C%20T%2C%20p%2C%20u%2C%20v%20%3D%20ucp.N%2C%20ucp.T%2C%20ucp.p%2C%20ucp.u%2C%20ucp.v%0A%20%20%20%20%20%20%20%20%23%20Check%20if%20problem%20is%20feasible%0A%20%20%20%20%20%20%20%20if%20ucp.M.getProblemStatus(SolutionType.Default)%20in%20%5BProblemStatus.PrimalFeasible%5D%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20For%20time-constrained%20optimization%20it%20may%20be%20wise%20to%20accept%20any%20feasible%20solution%0A%20%20%20%20%20%20%20%20%20%20%20%20ucp.M.acceptedSolutionStatus(AccSolutionStatus.Feasible)%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20%23%20Some%20statistics%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20print('Solution%20status%3A%20%7B0%7D'.format(ucp.M.getPrimalSolutionStatus()))%0A%20%20%20%20%20%20%20%20%20%20%20%20print('Relative%20optimiality%20gap%3A%20%7B%3A.2f%7D%25'.format(100*ucp.M.getSolverDoubleInfo(%22mioObjRelGap%22)))%0A%20%20%20%20%20%20%20%20%20%20%20%20print('Total%20solution%20time%3A%20%7B%3A.2f%7Ds'.format(ucp.M.getSolverDoubleInfo(%22optimizerTime%22)))%0A%0A%20%20%20%20%20%20%20%20%20%20%20%20return%20p.level().reshape(%5BN%2CT%2B1%5D)%2C%20u.level().reshape(%5BN%2CT%2B1%5D)%0A%20%20%20%20%20%20%20%20else%3A%0A%20%20%20%20%20%20%20%20%20%20%20%20raise%20ValueError(%22No%20solution%22)%0A%0A%20%20%20%20pVal%2C%20uVal%20%3D%20results(ucp)%0A%20%20%20%20return%20pVal%2C%20results%2C%20uVal%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22The%20solution%20status%20is%20reported%20as%20optimal%20because%20we%20previously%20requested%20that%20solutions%20with%20relative%20optimiality%20gap%20less%20than%20%240.6%5C%25%24%20are%20treated%20as%20such.%20We%20can%20also%20plot%20interesting%20statistics.%20For%20example%3A%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(cm%2C%20demand%2C%20np%2C%20pVal%2C%20plt%2C%20pmax%2C%20repeat%2C%20uVal%2C%20ucp)%3A%0A%20%20%20%20%23%20Display%20some%20statistics%0A%20%20%20%20def%20displayProduction(ucp%2C%20pVal%2C%20uVal%2C%20pmax)%3A%0A%20%20%20%20%20%20%20%20N%2C%20T%2C%20p%2C%20u%2C%20v%20%3D%20ucp.N%2C%20ucp.T%2C%20ucp.p%2C%20ucp.u%2C%20ucp.v%20%20%20%20%0A%20%20%20%20%20%20%20%20f%2C%20axarr%20%3D%20plt.subplots(5%2C%20sharex%3DTrue%2C%20figsize%3D%5B10%2C10%5D)%0A%20%20%20%20%20%20%20%20%23%20Production%20relative%20to%20global%20maximum%0A%20%20%20%20%20%20%20%20axarr%5B0%5D.imshow(pVal%2Fmax(pmax)%2C%20extent%3D(0%2CT%2B1%2C0%2CN)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20interpolation%3D'nearest'%2C%20cmap%3Dcm.YlOrRd%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20vmin%3D0%2C%20vmax%3D1%2C%20aspect%3D'auto')%0A%20%20%20%20%20%20%20%20%23%20Production%20relative%20to%20maximum%20of%20each%20unit%0A%20%20%20%20%20%20%20%20axarr%5B1%5D.imshow(pVal%2Frepeat(pmax%2C%20T%2B1)%2C%20extent%3D(0%2CT%2B1%2C0%2CN)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20interpolation%3D'nearest'%2C%20cmap%3Dcm.YlOrRd%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20vmin%3D0%2C%20vmax%3D1%2C%20aspect%3D'auto')%0A%20%20%20%20%20%20%20%20%23%20On%2Foff%20status%0A%20%20%20%20%20%20%20%20axarr%5B2%5D.imshow(uVal%2C%20extent%3D(0%2CT%2B1%2C0%2CN)%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20interpolation%3D'nearest'%2C%20cmap%3Dcm.YlOrRd%2C%0A%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20%20vmin%3D0%2C%20vmax%3D1%2C%20aspect%3D'auto')%20%20%20%20%0A%20%20%20%20%20%20%20%20%23%20Number%20of%20units%20in%20operation%0A%20%20%20%20%20%20%20%20axarr%5B3%5D.plot(np.sum(uVal%2C%20axis%3D0))%0A%20%20%20%20%20%20%20%20%23%20Demand%20coverage%20and%20spinning%20reserve%0A%20%20%20%20%20%20%20%20axarr%5B4%5D.plot(demand%2C%20'r'%2C%20np.sum(repeat(pmax%2CT)*uVal%5B%3A%2C1%3A%5D%2C%20axis%3D0)%2C%20'g')%0A%0A%20%20%20%20%20%20%20%20plt.show()%0A%0A%20%20%20%20displayProduction(ucp%2C%20pVal%2C%20uVal%2C%20pmax)%0A%20%20%20%20return%20(displayProduction%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%20The%20horizontal%20axis%20is%20time.%20In%20the%20first%20three%20plots%20the%20history%20of%20each%20unit%20is%20represented%20by%20a%20single%20row.%20From%20top%20to%20bottom%20the%20five%20plots%20are%3A%0A%0A%20%20%20%201.%20Production%20level%20of%20each%20unit%20relative%20to%20the%20maximal%20possible%20production%20of%20any%20unit%20(%24%5Cfrac%7Bp_%7Bi%2Ct%7D%7D%7B%5Cmax_i(p%5E%7B%5Ctextrm%7Bmax%7D%7D_i)%7D%24).%0A%20%20%20%202.%20Production%20level%20of%20each%20unit%20relative%20to%20its%20own%20maximal%20possible%20production%20(%24%5Cfrac%7Bp_%7Bi%2Ct%7D%7D%7Bp%5E%7B%5Ctextrm%7Bmax%7D%7D_i%7D%24).%0A%20%20%20%203.%20Which%20units%20are%20active%2Finactive%20(%24u_%7Bi%2Ct%7D%24).%0A%20%20%20%204.%20The%20number%20of%20active%20units.%0A%20%20%20%205.%20Red%20-%20the%20demand%20at%20time%20periods.%20Green%20-%20the%20maximal%20total%20power%20of%20all%20units%20active%20at%20any%20given%20period.%20The%20difference%20is%20the%20*spinning%20reserve*%20of%20the%20system.%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%20%23%23%20Reoptimization%0A%0A%20%20%20%20We%20can%20easily%20add%20new%20constraints%20to%20an%20existing%20Fusion%20model.%20For%20example%2C%20we%20could%20generate%20a%20list%20of%20random%20failures%20by%20declaring%20certain%20%24u_%7Bi%2Ct%7D%24%20to%20be%20%240%24%3A%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%23%20Random%20switch-off%20of%20some%20generators%0A%20%20%20%20np.random.seed(0)%0A%20%20%20%20def%20consSOff(ucp%2C%20num)%3A%0A%20%20%20%20%20%20%20%20N%2C%20T%2C%20p%2C%20u%2C%20v%20%3D%20ucp.N%2C%20ucp.T%2C%20ucp.p%2C%20ucp.u%2C%20ucp.v%20%20%20%0A%20%20%20%20%20%20%20%20ucp.M.constraint(u%5Blist(zip(np.random.randint(0%2CN%2Csize%3Dnum)%2C%20np.random.randint(5%2CT%2Csize%3Dnum)))%5D%20%3D%3D%200)%0A%20%20%20%20return%20(consSOff%2C)%0A%0A%0A%40app.cell(hide_code%3DTrue)%0Adef%20_(mo)%3A%0A%20%20%20%20mo.md(r%22%22%22We%20can%20now%20reoptimize%20the%20old%20model%20with%20the%20extra%20new%20constraints%3A%22%22%22)%0A%20%20%20%20return%0A%0A%0A%40app.cell%0Adef%20_(consSOff%2C%20displayProduction%2C%20pmax%2C%20results%2C%20ucp)%3A%0A%20%20%20%20consSOff(ucp%2C%2010)%0A%20%20%20%20ucp.M.setSolverParam('mioTolRelGap'%2C%200.02)%0A%20%20%20%20ucp.M.solve()%0A%20%20%20%20pVal_1%2C%20uVal_1%20%3D%20results(ucp)%0A%20%20%20%20displayProduction(ucp%2C%20pVal_1%2C%20uVal_1%2C%20pmax)%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
116bab5d00c38a7c22a21f9429b5680bd87dd9dfe8c08f480fde75ed30a34d77