From f984bc4cf84ee7120c85f0672001d8ee7c6a3cba Mon Sep 17 00:00:00 2001 From: Nicolai Date: Tue, 17 Mar 2026 11:52:27 +0100 Subject: [PATCH] Initial clean commit --- .gitignore | 14 + README.md | 2 + src/optimization/__init__.py | 3 + .../__pycache__/__init__.cpython-313.pyc | Bin 0 -> 287 bytes .../__pycache__/model_builder.cpython-313.pyc | Bin 0 -> 108785 bytes .../run_optimization.cpython-313.pyc | Bin 0 -> 34471 bytes src/optimization/model_builder.py | 1479 +++++++++++++++++ src/optimization/run_optimization.py | 731 ++++++++ src/preprocessing/exploration_preprocess.py | 505 ++++++ 9 files changed, 2734 insertions(+) create mode 100644 .gitignore create mode 100644 README.md create mode 100644 src/optimization/__init__.py create mode 100644 src/optimization/__pycache__/__init__.cpython-313.pyc create mode 100644 src/optimization/__pycache__/model_builder.cpython-313.pyc create mode 100644 src/optimization/__pycache__/run_optimization.cpython-313.pyc create mode 100644 src/optimization/model_builder.py create mode 100644 src/optimization/run_optimization.py create mode 100644 src/preprocessing/exploration_preprocess.py diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..f3b872a --- /dev/null +++ b/.gitignore @@ -0,0 +1,14 @@ +data/ +*.parquet +*.xlsx +*.pdf +.venv/ +node_modules/ +.DS_Store +latex/*.pdf +*.log +data/ +*.parquet +*.xlsx +latex/ +latex/*.pdf diff --git a/README.md b/README.md new file mode 100644 index 0000000..089a491 --- /dev/null +++ b/README.md @@ -0,0 +1,2 @@ +# LEAG-COALLOG + diff --git a/src/optimization/__init__.py b/src/optimization/__init__.py new file mode 100644 index 0000000..2c68ca6 --- /dev/null +++ b/src/optimization/__init__.py @@ -0,0 +1,3 @@ +from .model_builder import build_model, load_tables, solve_model + +__all__ = ["build_model", "load_tables", "solve_model"] diff --git a/src/optimization/__pycache__/__init__.cpython-313.pyc b/src/optimization/__pycache__/__init__.cpython-313.pyc new file mode 100644 index 0000000000000000000000000000000000000000..d65bd4bd47cc950422e6fa03dd987ef5a4e7f9ce GIT binary patch literal 287 zcmXv|Jxc>Y5Z%2^h~_x(C#1c`#mYt`B0)GsBf)yuyIW*+?>5YC1k?E=#NXi0vE4?9 zm7Q>nrE^hV@n+`Dn~yHb#*Fat4Q5U8qc8rE{3YFll4rJH3x3WDj=V@PDFliMlcZg9 z9NZ?k_~c~EMNG?;u{AWV);6ZrsTzVRt@bhX&OYeh&SgG~+13MFzeW@^kAt3uL~#V5 ztOa;w9hEH34EfENw*B9F>4i!YbAhb oO#s+mKE|CRy1KIZ7*pEX|LpgOz69s|olW1^bTdAD8sRR9Kl*`B4FCWD literal 0 HcmV?d00001 diff --git a/src/optimization/__pycache__/model_builder.cpython-313.pyc b/src/optimization/__pycache__/model_builder.cpython-313.pyc new file mode 100644 index 0000000000000000000000000000000000000000..4a9a06448e92590b313f56a379d7b6ed6e61e8ce GIT binary patch literal 108785 zcmeFa31Ae*l|MX(E~Jr02a-4pNFX7^EhO$!U~U9RJ&k09B}UR91kwn4Bo4&~i5(k~ zY!FEtgmYMSoLG=-ER$?3h!caI#Bz2WPg@|7UMK4un~k$M@RI#@oz3_EzTd0v>h2l! zNW#Y0&VL)!T~l3M$9u0{z4z)>)s4i&cpI)O|L02$_y4EO_AC0quNXC8zY%4(*&eY8 zHo?Bxw#hEpvncG??2sIIIyXBvMM+UwT(lI!-o;9B>^W*vyp+J+MQt9n$tAg1IC^v9 zrX(o|;TXv+j22?Y+4el!F&!*887 zTdh_kY;>_>Va2_5dxFy)KY@$8KmhafwHfs9p177K1tFI+LT8RlXOD(=SC91x$ zxkhmus8d{$&s$fmrBI>*Exx9FrxI7^ZShvuHAqTKi?^oH7f@VUlIlGT_(`qG%F0r* z15)jr<^wGaO$}`fQgfP`>wJyXHLVSebv|kKful-7V>1d+bCFW;9aS>06t`oxYg4Bk zTle~yi5*8CyZ6k!U6p;47Wa=?a@_TL&eW&2J+ZBOOJ7d$@t6U3mh8^zi_3n?X&bxr zjf9kAB{%T7t=2(z?b7k^J3#*tv|g|;sj%6+lWaB-q)xBPcv{ob5Bd6GRoKu`Br^&-nT;(VZEYZflt7Hs^;)PL*9exL>+{IHF zPm#gg?-ZS7Nk%LWP16Eix6M{bRKYGGUbH9KgoHvJ13o*xPv!58r{E%tg2EHGD)vpv zsQr?6Ps?GSv|n*-DpwrBcEz!}P;ssspQ;@U>3Hr<;8ZeqYw==2m}Q zAZw=&M4h$g#pk4~%{-&5wozI2zCC#F^Gif&`SFS~&>|^uwaxw(Z-YOeIQ$2c=r*6! z98jE%KED!G)7;#MU~{b!?Q1&Fa#V@mLEL1mBsEJ)R0Ce*M=6fFJxcVx>ixc>N>r0K zuwQXX&4&Zj^U zA@<1@**)`k)En`xlT%MjZ6AGN=CL)`V&fk^`tZ@$lg0$cFYod7lwaBy9J{qIX8#G{J=NRtbSeY!w*32Ka@9^|LerVx*%| ziBrl!Iw_ZKQ*q08O2qNN2fa1mjJ_k@+7@X#VmK33F`lGh{E$y#xSoR)t0>hn+<@zN zSG!A2oX|O5PMj&n&k6E9nBv43d^JQ&sZ90J8tDBG1 z_@w%*a&U~S!m!^(O3Zp+z}tiYQyPt|q)fVvq1!mz6ld|4ElOOQuc1Zqwl+!Wh-;;{ zW|lI`#8TQES<@7IzDgd^(gezABHh9WC591-MDUOn;--vJNv6pQHp7A{-H>QNtb#L( z6peIh7&BwkL&NKO`~~jC?U?P2#PQdXCTjoPnQtY;q{JNCcsljPFS?Syx0Mmg4Z9!!*~*BN!mhyr z6{o2>09%YI)2Yh%&^)lIY=u)8e`eXP9b#~9zsz>1pM3ZI*Cy{H!tumM)O~F3E(gIQ!l~9 zQ7RL8tiiSXF_jQa)JxRU`eTjQirVQA$Qe8|TE#bdP8ET$cd-@qMGgooFI1bcN*qN{sBMg=xzQ(i{cTzPhT6R?S#`cW4gLn7)RG@3k@N*&BH;Hb z&U#;q65~B^z~`@1qW3`3ZeatD;%xFB!L7lM^fcB$@~QSWH!=DC04uT@lKz2!64TfW zsR_?kWP)e3w;EruhxbSz4xchcqw|`?>Xzol>OE3(6XP@6kYES?0-5OH$86VKw(%2= zt-Chv4v6&$i6<9+VqtsX>4m2jb}Z~mn=+7)BPZlsi%UA0cp~w0zK%7Wv1hjQW#nIp zn;x2Y!Dm+vButSLrkK;loY~x$G3`oRerVdIeQ9|E3De|+Y39TU@*Ok#GG=v0%Ng^p z#4QN@WJzDz)PaONIU&#dN&J~@eHk;mTIGzQD{=FGo{;k83=&$#7WO3-1>=f-?Q*xr zd~)-#wXd54(d})WcgmU5`ds>%9`0>dueZ3HKgZb)ift~Wm6si|xbZ`9T| ziPq0TU#{~W4ZK;0ANesVNsjU~H2I+XwRn|;hCp+z7b0<;S5lk}{uU+bs1L!Y!#>}B zsF%@DFI$_G*udU~J*190kv>!$${;dq@Yh$<7b?Bs8=rJH3OsjHyY5KAC*WyzK;wQ-Q6Z+yO9b5ByLelZT$s;F@d~$h5p`0-B$+S-4sfx!d zzA`~hm>rCoJy__ZUeA#8=_58V_gK_W;WWB$&hQd+2mu-U{4 zaT+m`a%e%KWjg0YHMC1_R}{u`(viZhvb3-W{H!GDja{6Kg;7NG$RUHjF`gsKfjR~;ftSS7h9v5p77YF)Fr(cZ)(Ds;8yFbD71s}dmUqKE`5XR173{e=j*zWmVV>*9bNeE&W zVToa44v#UOyrpFYMjQ|EcN~U(*RZv5X|>RCKHL{*=lD{)_7y1=>>9yWcj#Z~bTfjl zCg@+04$CfK;v9R^RbdkDS4=s81BG9cZJytikq{AfA$DF^6n-MMXmG5MC1jgQA)Pkt z5+*N>0UmXVHl>bKr<*EWn6fwu`DgKB-=qBXo;HH~Zvjlvs%W%o@=dhr)|IOB4IxL! z6{ZS#!nD~@{zPkPNwe)p3OMH40uC6v!F`+mD2!cl*k-GPJW+yCMNF#P$XikI&H0){ zUF|4=W=-a7zuTV#4wTP0P_l@B;XKHnVv3#4Ig>xt6gxwYO{?6@OA}LsnY=vz=#cnX zdVG3F{A@iw!xW!t$|o}<-pFT6NW78HSW|qODW7p6@kT!5L*k8m?tsu&!dqrMIZsQo zt#@v~=-LYS`X~4&R_c6HoTze3&h?s>BPNZY-L%59+`6wX#ar)1^!XLmy`H^42i#69 z+h){-hlDwm+j+Rm!mOmjc#=t?@bE8;Cq{!wyw-fMr~^-}Eaz#A_o(M89&0?&z7<6d zHXd8obnXainl)@qvxcdu{vI`5ZPbv5Of^j`yNkyfPZORU)XNy2vzwAUNyyKl=Q9jt z0fR$f7}C~I7-l_qaU?8@*&~FVc`Cvh%ag-l=Ox~8C0gOc!mLw<24qc9 z(=E#zu4VdH;pK_Byxz?6G?o9V8x0yYbh`PTB-peUc$VY7(6a(l68hwx&t6FyF1y5B zUOWGc%4$7!svbL2oB;~IlfN;ZArw9nY1SFQ)k5Tkg%9ZpoRMH zCF5-VS*CCH@;AowcDF#BB`jsNsBsCN(%nS|kb$I3d9#@sj;>O_<^x&p;G8Pt4WyEuAvNxgk<1 z8Um>_3o>V{T7tg-GUsq<=8YgtA=02#8!$=|Je!grc^T5lL}8^FhGKyyiN>ao6stl~ z%(F;Q9Fn3?O;LgYp7myff5enEo4}$4IXUUjEE0>#4h~KytkzS^7v~FWNPh4yU~S8? zZkxD(BZEe3VdXZvaA!zOwbqtQ)@o|l3L@!!jf3kY7K}i@B>fKzTgWQCJ{wS_Fz$t4csqymZ864b ztYun@K63h(iObkXJ02tLQmc_xT&C(doMU6Gg|(ALgA&M$^u%``53hq^8K(CbA1CZs zzGFoIv-8YygV~8z%#L8>UG85YE*GIUZ+Da=s`>d>B6fm*1)LyPRX)sf7FY5bpndRr zU=l?!8wZ!8t;L$7o&!=t3fi)=yja*c3u8Qu66LGZR?~P`k~L_AIL^y3o`y9Uw6K=P zo2X#8Z**B%T&a&-sEW8!^%B7NW%I8VSJ#|H?N+l^{J!buz^6s8M|(dnEVR@->Y1?)0yX z1kDYfL^*HZ?I#S_@sJh1kn#wBZ#+G3XxbI?Q;3|V{Lb>VMJ$a8rYpoZ80ME3MGGYw zRCS5fC}S|@!*FBHXX{PnGZ;I6r(&n%lJh#R%KZs9qZ@{c4nTeRh7%CyFpsa$LM znB=Vzi&_ph#%|@`M8d*`-)iHg1=KJGT*ceSGSmxrOEVG{6o_x1hV-%tZ@ag)luk=tkZ8{;`B4P=4vi(BDH))B;6kgiDbS@I3@*m%q9I-+`vq)mhC zTJnrGLi|8`4rtHcsZY#KJr)iFv^<_m0-EQEZmlrVatno}zeuN%j)%&!dB}M3JevMs zqU>qLdmfr>GrR~4PYBN+Mq)Y9TcY|Ul=3_QPm4Z!GyEpt@Fz806YT!YrrPUv23`&= z6S6J~u@=3FX6CiI@iLyrJI?XXf z{*fW=1kbOyiJC@kGv##GP2@Cmop#=`oGNc3r=eT7>z3uD&k6}UAaQIlVenl`aN-uo z8-$B0UR&`GMRuW_bOwLPJ1)rrGGd8VlFr{pkh=xe#cQZv1=SDv(R|ZxUh@&;XH_rC z?=Dk*UY?)s>1OzoY}MloZ<1}EYmh8RcO+fkqQ-a&J+W!INzWtgetD5o*s1AeMZ8So zIatT@Zx^?tgeQ$OJVcU&>cdd^LmpSY)srMlpAUNtU6^x#7Oh_5L~%P75mqjCGL3kK zA9IQ2n`t~p?p-{U@icN5XNsFi(?q&+#(N$z(v`GnBfq9QG6-GJAY_39K;suSj#c z-wb~boaiN&0Ef)d_g(dArU~J06HRzbZ?@r$QAZJFSW_}MZ3rb39Xo_7M#ty}gCMm- z+=0BhE$X|R?0N=OmV3k<;jqxD^4pQVH&NtwBBdz-Cnk8H9@ez?*$))$y(#_<5Um~t zqUSAI`(_wK#la8-407ef7AI{V9r#(##@|U6VJ~Y z%NPkUP7h(<04HSjja57!B?IIU z86V z{V?*ve}16kh2KGcX&`#BlvG1#3n}+SOleI;Wz@^-NWRozc2YOl@i3OfZTq9oWIUA?paGdH!*sPb(-|z z=dcp5prl4tlJAz5or;Gpq?i| zs+r}v_Z{ZBSIbi}@?_UkzNb@Y9!v{Ra~8l0t0=4BaA z)dyPmL+KQ_h&%Q?If}{7kp)@|29Hk62PPm0-hd69u83x^`$zLGPAeO{dz6x(S>Z|B)8AjebWZmD8b~I06JUvO}A1i-A?{T{I$ecp!jh6jY{FtbH zZQ<`JYW?vW>a*gA@afqO|AXR#!{Fd?4FVtcKO{Z`&h!Xi^|1KxaF8wkhz{A~;&BsX zCo{+%9}#3L8Dx_@{W@OYcl&{dG(5oAY}9G(`4#nP>S31pu5fY$D@n=wOJ8{;w1Zom zFy(7$dl65d-_g!npW$>QoWggX816gJJ2-pzh1ijEnNMgv9zGlD_1>2$d}f?YEA?}{ z)b;kt89eN{#>c9r0aIKf>FfHou-f68-8!`{G3*qHF?CIm>2(dIW z@L%HJY3+H8rxZRv&errt;WX}lGS$M8uiQ-E5dVnQbx`Y2BJr?JUg9k+0c54h47W8e-oZ;`lV?1w_|?stJq)ml&LM#`1i&WHm6hA zrTA-XZNcC3c+)qxna(iU@8q>Go~HUlrnGPJ^rqAwp4vZcQTyXO&Ul7c&boQLsl~&6 zgH?;_g+6W3$FOUreFJ6y+0FnK9E`mOAB=2U!}p}mM_jXt_bc#zMfoX~SD~gIf5!h= znoA;!nfRF^aJ{m9yr#yJgflk3KM`82mb(12!={F?$aLP_Mov5w(YmmUr!by+iJw84 zZ0Bq1TJm(%4$azN(H1=w>QAW#TO4j>oyE_HBdfDkn~}9FLY*yoT|_un8Sv#H9Wtgm zgD(t$5u2f{;T*K@!!PX z7R3usG40}WMsGu{bkBZXdvKF^>Cf9wVH^xcG>Xa|CVcCZsQ5&|=IM zAJH_Hf6z-lhmt>k8%qAXR`Pdv$?%2xI#Is=w1|J&9>>N?vkuK~>Ko)YWvZX0T!+?L za5LfWnO6WwpwH^1b`0aUv_a!Y|M`SH59gzRs^^hWQ$Wc?^M${Jeg~ z)AI&-otJFVXp&YcKFuibj0qax#o^ZwbL>xqr-hC84Ddn_>f!H9xqqI=8PA)`UFgVFWvGvA?p+jLor$<$+>XZyszu6?1f_s^QX_-8G~_~NJHCGoQO0{jPi7+$_G zf|glWUO+n*_L>O&_X}Dt`mPRx=TQr)!v$7{FW$;Jh+oub<2k)FMf_*+b(D6ImG-4u zSz3e^e@Scc_w+jaR8&|UzRT+H*G5#**TPF;v#ix9?Q5(w-9N~> zw2~HGN`=fx(yZ_!NN!J=dTFH5`dxAs@hOc1|A^O2cnKwU-G-98w31)qB}ec9>gBNY zpB2wy59m(;s%M8fsi)~*gg=V^0*T>e8bcs4oD)?#Em=7-u>Mw*}&^O~MU zW-)!H`uvj?W0dueEar-ERsZS-ET-QS^P>0?pf$io+^^pT%zs_Oe6NoAUX*{0mH&;~ zQ2sZx@?X}=|1rw{DJ%c)Z$tTiua*B}z5E}ce1(<&&D&7^H?{KrQ7``|DF2^X`A^@5 z@}Jhq|0iC)smvN4XFQo!@-F$DsRRqHUYj3%g_kM34j$NT`sTgjzTFxme!@!@Uj-(9 z#xU{BZ7Ba4t^8N@@?S&wKWF8izYXP|*UAs-<-daR|Am$Rt=mxkx3uzQz5JI^{u`|P zZ{LRUzpa)3nqK}7QT{Jj`TuYm%Krzg{64+>AjPW!(@&{uBR~V zVLCLz$xM6D=Ylv)Fy~>?RZ8HG$$aXdw5V#)BQ94%ku~P z`nmq=cRASm9_P(q|IIyWh@Cv6g^gINP?UcK?TJ4`TBpC%zg*8ZQuto~J$r2Vf7&Xc zquAQA>G1CAKtr1k$8r1pzCcz}YoH~o#+P-VIndD3aLA`5NzHrrH}7rq)%b7-^Pc=z zX(zsvGU%2`x3P2^kDHRR!zaYMX~hB9DiqpAQJc$wAfFm6Y4&?t{^HhAHfu8>`!7ixVe0R;r1idlx2Yv-+SY9| z%6H~VyQv~D_Sqzics2-TQP2}S-5TXHTZp5^)Eix*wRq%^Jnps=n>rW;&W0~ zTii^w9$9TMGYJb>>j-Ie%}w40f34Zl7K_|qH?rSYg#aCMP6QKvA(luw?5aquC=#Q! zIbe`Z4vrk(n(t7O!wFxp@6V4>qMMpQzlyy|PO!TfvSFN7JKNrZF|G{{21)2raIUuc27F%Lg z#SA!Q3m{9#GLWu+Ivq`9c>nJ=f%;&^;`2hcz5U?nBd3mZ+B*w6y$j0&a)kqg7A0W0Y^;s*8DIDDt}*$I}A z;;&~DfXH4?y*Byurc;|fdH=P{F{kf2bx&t;S9Z6(vqa9E7fhaaEp7Da!c&Fq2RlkT z3)=6M({h8ZTormmnEcW?7eS48r>ZE64i$wi75#3z6cxw1)D{F*M4P%42I(QYRE)>r z(b%NcsW764btlA3y4sPETMn2I(Onv%9#x_4LtGM>|tGi#rc?Wp|Z!PUz2^b1_rvRM|Op>C5=y zY_wW0>i6L^LU^K1rS&S8xfV+HuoWeZAwXi!U$wqe9Ly;0Nq;``VrE}PabH?-&{eDw zgqhs+?(_lWuPS#H?Zfa_kV32mFZ4x54h;1=lU2cP#b#APm>mJW5x6QXurhPi)fe-+ zOF;tnp1QZQw9DRA(BjG<0{2j(l_i;gK@k}^ zz@jY3AXEf2W8Pp-)Jsg5XFN*9eAsBkFz=&w=jiR8uXb<&1Okt!VCvHD`Ohvq zzpyWLX`g#(FmCDl+H+~VU_JL{Ogkf2J(9yAC_PJEfJ#+WPD6TA77hp22pY$IVHjV3 z@o-Oe`}~2ld^s(j6i%4&3c9BEr_H|zmPX_{fSIWutC)T$YcVxG(mtA*sv_5wKnfoCwm3AGLCoPm+ z3$M9S23&X0D;R+~%4JtBVutzVdTdm3{B>K@sCb4ysS#z>I!*AQ;hG_l`i^Qc3NaU{ z)>r4dH5E}Ytcz|yV1;vvO}D;67K4uFw8ATlL{SjD*`7<9{soO~keNuL=`fgH*pT@! z!r|gJ70!C9siPKJqUmTg%@jp4ww@0y&5+oAN3X+W{xdDaiqY$d)hg6P-A+2%sb6`J??#e^e!5gu zQPw&HRya(Qg&C^Ilm)tSPz^M?f-Dn8S7$Lfq@n~C9Y$IYA-O?XOf({@G7(iIdt+0x zzh!SEs_G$BPDI`sRV_K!)&)JJqZ?X$O##JT&S(qL zA~Py{oa)Ey`BYQI(bS~G6FCQ}!QGmaL`K>z&6vbNZJ6}9dp9vOwFL?+l14ovS|p7m z-=TVMMoA#SYq|3Vau>~A@xsqP~TJDlFYdea2YA)1&*4vxi+xqI^mk+-@^R!pa zU2;5LPA><^UQJzcwyiI9NuPU3FmB1)H=MSN3VXnf{7xqnXU2Rj3AdgY+~^c+{UU5L zA}&IJ6LFbxr~M8xL^`cqJbN;T{1U#s2|+*x_);7dm?9|>+R4GNy%cF&LQL-gCypO6l7RK$od z<7<#e>?U=VTH{0-A;j@w73Tl{NLMtni`4n`&xi zR6ZDc|AEF@U{87*@!unaV`ZnypvUa};Y`qI0E%CMYY_Mem>@QIQ z)AYlJRc6Eb0YM@QDO5qj)Vs)o0RvD-OBhR{ z=0s}A8(}S3d_J?QxP9&El2av}u!MK!^pDOvpQ)n6TIzWcxkpw~)Otr!QPj6XDJc5) zC?))=NY&rQguk=d_}g&RBnmB7;#JjzMsDf%_?AULDPa-6K?LI#7{gT=GEE!eLllx} zmP3n0%HI%-cOp|&8QF4TOGiq_XxTNfGkqYtK+Z0JM&`#(Ppjj9-IOhbbX{ z8~BMO{3NLOVFWGNcnlDvlDJvoClU|E)@YUp{v9Fc%C|(9RnqIB_$oP3(vjVfBfBPb z<_=6QlqVN<7kA$&PhKp$788by_o2TArw+%Lc_62;g$ZEPo=E*kv{N$c-lV>Z+?Vuk zkwf^Uzq9B^M!qHmUva7IUUi~Tl8x@9#X#uQdIFqm*n6#asqdk^8BCVIzUss2J3+ktFY4;yezMCa-_D7&&d*AC2IaFDbY^3SyhPOZwAp%v@!R!o9!NFWoorCQRBPh{g z2n#Qf2C{tTR>l@WC^dAd(risdMI$-hnLJ*itSA+yklNtuI-JLgoxEP*>~atw*~(@s zxl7iAh*4egn2tysSar^$rb0%AufVXjvP0E}$>JW#R(2z+Iz!;xSY6J@arf@s+*)R?X#{JTIY}|q+>4t34~_kaC)b@=*Iv->Z7+h9UGUiQ zGt0Zu;K-8FHC7&9D7y+@ce(rBE6>*UxmWhl&@$`?VI#`_4c&r<%nw$#kmeCvdI9y4 z3K8O7C?g|$HUkxj@V64vBOg%a*8~csVL~ETRn_xwAwrgPyLB+#8U;n9TO-M4DwzBs zRWp+lZR3~DZAatv3J0G&8UgP^p5d3ljfvMH4;sJ2u5kWXNi|@X7|M~OJWk5x9r5Ln z)Q$LW*cDzk&btShhu(MIJt87qU@03Ps`fTE4#l4<06mlAiF@t2hOU%LDKDk{Fs-N9 zl0RQgBmUf=@#kp9pDl?yI*chxD~Tu<;l^E#wBbcb3azFxWdv;zsp2$jS+eH^XBMbd zrZdZAS6)~4)3eUb>WS_-cqzLlL7sl6?7EZj>$P3A-D}P@^tspexz-Nm*~8}&&cw~M zNP{R-y#Ip0Ff$*P>mn?BgH%#T20hliX{re(i3CLlLWoSFVaZD$I) zNzADjEcM84Z~46I35YIao*gU%j^ZsRVI^eZbJZ|H;7q2jSPFqmNW^5QAAYY^ELas1 zE2e~a!;aIGn%UZOwWi!KQIWA5%zZQEriEJa6mB6hBLnYodZo)Qb0?6tq@e8mYLw-8Z~34$*>!DzTNc!Ox|f#SI0gisaeJ*6gTBQdMN#KW3J zvz*dqBv2>oqSRrXfTp-qGsq4`O6(VtRLQo)XOA8t;4xlE1GFxrCXxioY*5S*X!7tx zT20WcMIM;b`?bq`)jhq_`&8`{kVB{UxiAazbk*6ao;8dW%4n5m}nwpeT!Z>HQ&e%?dF`QN@PWs$vl0Q)YJ3c<)6Or^jKf47@)Ai*9Ax2PlOM+tb*vJ&g6! z&sXYT1wHkKvNEwnSkcSr~NTjCWp=n97bp@2zrDKl$Fz1~-AU#aaS@y}$% zQbm_qh%UOKG>Ff1l5_XsKG))&;^((s*xEZkSSIu?l9zd8mxqZgWn6J9>vNT{nBhJf zkrcB74QA{QE3k)kv~=un4sx`dQc)1%Xz75Kga(FRA&!;~!_hL}grlXS$?R%r!+tgD zMha1u`XI6HuuGuDZ_HPc>4QXV?1DZR1z$T=AGC1WG}+6kR8hH@egQ`o z1}$n4oIu00K)nmRFH<+t93YzOicA#O6-J9m3?-6n<*>>B(t^Jp1Lkk!0rZ;Uv*L+; zohd!p&rkW@6f92G_$;Z67p9zYJpmsxnjRUfHpU|drZAJ&e#rO(4{*#K#4r_b6;lx{ zmqfDnq(Wk63@d8XRmdp3r#>BakqTKc1+-&cP#HmmEZH(hm1II$$O2mf;)V`+1BLPh zQejeCBf#=ef+gC^VY!=>$mf=HdwUPQ8h9npn;qPB_rR{|tGlX$z_onTs}2!hXW_n5pc{dy()!TiU?6OCtC)g;e@8LaMVd#OGC35DW`z^8dL+HZ4zoNvR3@wq>3#&~BnwfDJ13cp1~5V} z2PLSf^fd2qq1o^a3ql03LRl2S&PMDsNfx8(fJZaY>`q`i`MiKY+P zyhwYCVMypy%?K}J5^=`)om>WU6h~WF6aMBo;>XcpPXmI1Z3YDQkcr_3<9dpN;*J5a z;;L9dYlCWP2Wos*Yka}GtE`0E8XrVjS**bQ@;H)iF+Dqo!G6%mReBKMh&1X-1gj2W zf)7!uhv~*xa9o%?8_t4hjztw|NoU=SK?x$QX0b?CZ4m~ANANHWPxe{D*;B~wUCFg2 z&+P1Tm-M+xfZ&1c)$;c0;BH^Az9G1~LEgSkcI{)L?H+FF*wg3Q!^GQRzZ)dsYK!z% zPa_F;8cDd*V&Hp1Zdu0ukx~);JcL`mvn|EgmQ+$V_>}IU*Z1P4jAGlpVhQ$EQx5m5J-%#XI(kd(Qj)#}_DUy{_o#OyKzRNs^;XV_kNSY6CK zN1X5dUU-(O^oVyCV*5>SrI%^B52ID;o_HUpkTG#e?Iux)lcl=<);vp=%vQHYayhz4 zWSN89RA!u7crJi-sRi0TgxOexjtwT;@1@nL0c{hC3rGzOl2ztGsuru;LlmO=az_om z6bUzs#a1usx-#v_xzvIWlrsfEGR!mfi)2hIihQVaXsg~7N&btb|*&JmMhM6eYB zP6Q!y5g&^EMX_uuTHT4DgTP3wQc-O3r?kc}_A`j=P5C^*bTTs2FkE!riG_Lv=50ng zmh`7hJ@3@EocJ=*N{>?&2*|r}sq2nt+XStc2#7Qie4h%1hXLCZls1|u>b5`YJn!sF zE$VX@1>=fTESkHrj>RB8TR{%JF*S3bH)88H@=U9P#vNg!P`veygX=EJR4*$j$jb4A zZ5&E0KJ4NW6V|qnJP~#=vzF>GjQS&uD=Yzi(!N#dc85!*(+0uPF+@W2|0*7e^T5$& zC4Nn_KhPp!TjCZaO4#3U09#Qs9@xw5FiJFoq|||yfgeQbX-f2_>XMBW$|!9g)4Dx_ zjA2R=IBs<#+oYO#2UxpFg(iBDkass~`rk0YB#f35vM`y``W1hV{iSi2$Mr4RjBQ>| zuQ;`0AZ@CgHZ|y)ifv6#e&obQuxm|`?9S>ecxvI}3%d&7=ekgyQYgC%yAM8l_?g4- zq$=%clNWB1-J61On?B$a!0{J>28`8yR@f)Qr%oY*DM0!X5G#F|ZjaH8@T{ZXdqjGi zBA>u50zO+09KfdI!{F0M)ca2GdHBR(+K9HjO?FR)cgX(7_p^n8Qy0P*op&McQg*NX zQl7kIlRVXo&w5Y|t~XXHakP0fQJyl2h4wbo?+qEc_5Mx@@B7`~hr>31Ptr=H9*e(m+~&>VZvs6X>Xpqx=fy0##G&XvbZm}zcF}F3hr-}#Y6JcL$do&Fz(P`y7^GC zn~~EEot-`$?TnCaLRm{wL!C-CW-f%tA*d$BL^YIhZ$pGVyP+4bnuT)i zLTou2CFd@`WPd64a%``?*V`K_ui7r>ZkOHLgC23fvsd=)4L0lxHXjH!1Y}P{HuHfZ z7o8m5wKRHb$fbwm6DFjP&EUtu;L}T;FNTRX%|9J-_G2qUY-PxVlojJ0b9QgCj z=5=TH*t_%O=}TnylAhA%%P*8)D(x+}v{PO_5|U8THPQM|v$V&%SK8x)%#IUi(n5b@ zkgoUDTB@YK0Z~6<6R`OeSt!CUQfzLm>VkQ9I~=#!Nfs#6ieU78^dbv0l#5Pu7UL2t z!w_(pICIvqRU(Zq-+E~J6kY+Gh1D8Fi3eyxjr26Ypt!UbY|S_PEDWvZA$+Qu76`;y zxKxESV{mS4oSI>2Ww1`14p-_x0NFy-ai?LC-jXxetSa~Q%yFmhJ9Xbc=5#r8dN6tV z8_6lB<4(n4tBgB4vpXip$sn(R>BaK&Voea9zDZ8rbS-W4Kw2)&ZGgTsy0czRo86uA z?8Ij#UW&Tpy%aAOY?jkDzvZ;2Y_VT=Iun!BLFz+7HX$ToAVtb5g8t^<-e!5*0Xgje(Z@j( zeMp}{`8ZrC<}f?#N{Ov$K2qb8>Ki&?)DD>pB|>77@oEIBz;(>7ub9_TJ;E$~%LaRi z4L({&_tX3lvJaWL+QDY7ERWVv4kVB1n1V&lr9lj5 z1y2{Ag)wFJjZwD5X~&CCt~s%$Jx+U_Jhk7I_g0K8W9(b8c>9a=F>giVsicp2J3vFg z1B;5Y?axofP1~uGx&B~NF4iz~ADd?5OE6VE!v2(|6dR}|+0+tduKhDL9(7vQG+P#-7}8uL%YTZjhX7k^8PCGvSdoc%EB zX0r((fyrpUx+0Mw{Bh(E#!V9Hd`)x)78_lxckMK#5C|%ZkPEoR>^a?qU*CKD;K@TL z4q+X5-HCfTypQcUvj>}5dCoM-?&-m}4e(c8b79RTJI09E3tJE|cw8WYW?g#&0sxk! zFkx(f8v#Vdzz>?>Hc1EDSQSU4u_hXJUwxQ0cBiYLhUK${VKx4;3(GVCd+B!Ud^aFe z$^ujCYifhLYGtABwWW1Kw50KtwFIptVhXzqXe5nap;Pu#0*<;xChOFWy%j%uq{k7~P= z`%?bp{N7@6R$ec!_Q-BeFwQd=Zn{2lYv48<%w~b&%xa@k`YH92pP^MsbQ?JbT2P9T zOuH*cJ}=Gz)H;hb#KM!+6T}v45KBIJ^u$pJS_eBaNtiEvs`By5?v(C>-I-)cDt>=?t zfk={Q7$0>`wpGBRVm{RZS0w(YN^4|r(gt5Cn^>NCRS=W6{zNUJ0bFNohQx6X)prWs zfj_-G9wy&l<*vSscL#TS2X^m~ckc<-?+w=Pmv=Y5wz2U&tvc-!Fu3Y}9-`{^nyWqz zwnyyS+_!O8uzL4EwNI}01^3hk_w19a_rJDL+uT=!C&3g$lqT@)K_IfgH>5NvA;IOVxz_iaN6GiyuN0ZY;pBV>UQ#+;jGx zp5mT693#-)}VU8!)Dv}->X7qKWRK?68(b+aGw=&ToUZYm}tEaj2|mM zD00|Z zXze%S9W`%|XR3M{&nfD1R4vFsq25nyf9^-9*HJUn;WBZjCwr97OcG&5}o%UJ%Nm=Knyv>BI2gWY7caC35HiB&ArygoYYh|R* zQWZ#=!Dy=a#uyyyb5wzuc=^|3w$9>jIJ=6UUUPQM*W(`F@_KyYK;njA;szL|r<|M8 zm$;!Xe#3t_?KvCmD&doe&o62rBaZ_JMDG8_In?07R$8&?q-vzF)qt}>qiMZ%h4v%jdV%BvlJ#(ZsMTB zQ0Uw^c?M16#9fZ-P3b+@J5FBZk!OUcN0#u?tO*em1K#Yw2A1D;KA(2P~d57 z{t|(wI1+eDpBByvqzd(vc2AIJ(X8mS%hP%bz*u{i$g6kAvvve$l*{h&U|cyPsRr9G z{xiA}dD?8#`gb{P4YpCsZ4Ik7Z?QETyZ#qPz7y7<0A4cUM!1E-wBZ_$*TMu;CAzMm z$>$GY4}<|mrCZSQ74c+q^tcT6#g_c2;(VMn+wr7?rM_v#Mfw~Vkvl9Ob-k)3#aq_| z2Fm2HJowH=W;#P@bVbO+=}8RBBU;-X$)e{)$Ro;=6qZNKP3I9U#QI}AqiWWI=_Z8b z7N_M#15nr%p++&D_^>=i-E@s&w>eO!Sn|XrjEQXLHZh8Y*{TdNhJ|D4j=?BCbpI!+ zEE|U|A`oJG-pCo*MN>*t0=CVQK7{ti8 zkY3#?iCI+?3@BWe{5ZXOh;9$l?E%~@4bRrgi<@iq`D$Bm_FN=)-&csic#x?&KBfeI zJpYn6xU+g-r}yekFA3QD4-D*=uI`s0pSI_oo^@(gr@hnL87rqw>ze=c^0UzVdJgtv z$}@`lQ;X&Ol8pOKFYZ5{&lhHlDIOoXIF?|?u=g4QA=ul+uorpF%LazMce&;Tt70s% zrzt3?6H-`o2kr%G!Y)h~Y~4)NV5~-~2z?v4l>S1u|3kO`4>uq@Y)*k;l~{tI@_Yfm zq-3zsR|!MJaJAWisMMrij?HQ>zBYE^W0_ydY%jhZZA-}*NY0Xzvv9n_XgPUmSNcHy zVmV)RmdRf!C$GGeI{+t{)n#-diBK-D#v+f3VCC+C%KdUBbF^waAXgsjPd@mIv00cY zNRjNy*sQmi%kcxnImLOQ_;aB)41cT8nh^YbmfKxaSQ`)^9tAp-Io%6}V7qb^P5^afv*Il-Pd9Z1(4$fJ5 zY5q$KFE5nmY+`SYujTvO;~49C!AS-Eu6b|8+s4m)D*@lV6=%zs@D^sX3TftQ0?k|% z{v!Pjnz<@e((iaXkczrKP&{_^MEm8e6uQsP$2~*>(I+LpkF+7A37UjlV5xF2(gB0U zzz}(FgXxSkE%!EH-)At=jJ6_8!}(hy5SY#wfj}={&sEwJjF^6#3Bx=%kT|BKsARzB zA)tiwli5SbK%J;$oP~OeB?Y`e6rgke>o&*IcxCab?&99uSMy%UyHv_pM9tO88mz1j zHXa;k3|wstz{Wy&oFcn$`dwFc=cazw92s7kxL4NT{&F5}7c)8cA6Gmf6h+jB829)k z0wE}>8wN!*whx6O8rx}mib&jJ+DK4@?=3infkG4a2T3sWg^lecdj_LXTfSrr^xp*Q zl)jBy$exH{C>$Dlzb5RFH-grGV`dNAf3Kr=r4G1oTu2UjU_s|j*)_L&{y^atxo}G_ zogUV^Q!cFNcU9yAWL4;0TVo)e1SzrijT%9~aHwmZN<^$p?% z@z&pm_pAl+%97&l(%$*6E_`JnM66xjfnBv%chzF2kzkWF(A09ZsRf=N)Uyoi^t7Fd5<@Gc2{wCcGr42 zxv0B!V4)x{6jbYX5GYu+GajvflVr zgz{6FjH9Gr5NzYM`mPxn-{(qLLn$9y)UBNlo6zJTAi&N1n?( zK1rii{yWTnkw;@47OEsY2_P$@)zciSLtM&?*jE33EZtF;Zz%CvT70aBQ4=Lei>F$l zI*LpC%^IZV71BHNp8IViax93ywq)h= zc^C7L(*nw=kaC)goF*ctapR9CscYS?O(^&_ss}c6Z5eiwySNqU7ojPR`AJAWG3of$ z>oK<6#m^P~*@oVrHy7a>w&-O@6P=rc!=A1u*>bTqE*psxCtgd+B0y!oku>G4Xh&iW zBHic#&ctjK?aqB`R8(T_^>~|m!dnT^i4!1DyK~-3j6phQoaQ?Pi`*yC1(fIjT#zI9 z8pDq+qODl43l3(mVEb_iPIb;FKWZymz0D@I*$AVqK^`wK7pWF82!Ed6gIB~O6(7MO(mRg-YwM4){cn*l8jrpLRs~(}!wgXxfy)EnL*$=`83<38qb#-7~rl z%I>+rxVg+`tGM9nLrd1QDBZ=HcBL0(D7sZ-7TBbOSfT zO9(9AZP3lo+H7Fww@5KZWl1~aCgey*pTO zm$7L;r=4vWAm6Obg4 zl8WMN;@i#r89&vX=5}LzTctKrahV$v;lLAflBv`r5?F7n@y~y7Wr=j@PCx zJ+ltSS_dv1#j)0IO`F1b7Ctr_t8Wv?GgreyEI*Gdkmvz1_!`h}0x z7Hvf!Qj4$###mFLe`*9mcOb}rIaYW4~_hi1>n^KxUb{I1>-dG$^?XJ@dg;kBFwl?^8z z@=1-}x>kRE0G!x*^Q^BSW7Q?b1Exg^jA)~Q6~9zhWjp<%9T(}H-yvFV`#cTr=<=pu zK3`|nf>X35*#yk$C)o<4_?*AbZeo`XE8~-(4)DxA$-e>^W49`!^cG1AQCut7kZrTC z!SBo3VQy;Ho)@2!vf7+kbFx&4;M1r-I%^J+OSyDgOSh@GUH<%cHwlzsfrHosmji=a{7!u_sn42OhyW5J5xcx-Bdcz;;}-S2?PjA*(c^ z)O~8uPolLk5c{;wM|!$AGCk3{FKLJ=G!AJBX=r}}fe>E}l-QhNZ3@*b<7lA=N8G{Y z4hQe%$<~I^kkpoTw~ACxEct|lZ2OA-$HuvJgT}YBPN@^MRpI*-LWAtwifi?@iq&h& z)>q%TR`94AD;w{Q({~Ty#s@t^%TYJOSz#%lX{%?dRpp5Yi$qon550{@(!rMkmjjol^`+%r@$81-p<&-ZL!*rUb&ZJAhd^W39N9yI;JN$XW;=#HkXG!r zKT5Y}CgP@^4XFoV|sZD+ECJft>eT@SHck&tsevDe zP57nLZ#t+Qn^Yq`n+CuhGnTSlNP0`@j^zRKKFnn?R@|8NhiEK2G0T6@Iaq&Kr9CkV zO217ev$IjAe$j$vuqmNIokO+x6Q!6shPh8OwX+CI@lAfdWoYBR(_Ak)wTk~*JdrpkuaLoW`4Z_^g9u$oC||o}=eF{j9t4(K61=9K4?70z zA{?8>fc@w0!x-RwsrGX10MzE7w6>8nK%PQJ<<}MBTNh`RxKt1pLV?_%>=k?y#mP>S+fpj0;=r9WFc0L%|&NUYo4%Nf5a^JI}KG>8P(vg^Bt)U_*(d$un!^+)%)vV z3O!1<<8)&(e=I_<>_Ngci~QCHeUX^2al{Dv@8^grcOXkm+K(Qc`pXG<9rkMzCO#&7 zNkCxC*bXl|5R2QP7(YGb>=ZgL_v>@y%*8iiZRvTcN`QqGEKNsZ|LCc|m_W{k>3Pb8 zyti)xRg!jP2kVZu~eljFY=b+xV?-b6%>MZU&*qJG3`yYUYPVM$@<1Fnxk0eP-iLV88T0iz*8Yd7Y`fH z$nzA@0l3FtQ)JpUwrvdPACMjZ-7;m8EC^p9> zO%}cqEmFe)%ZZCvC4psA?K|SFZBgRc7PTzQ22{n~CczGFW3v|$Mt8mQ(g4d0YDyAO z>K>`N3En(4jlKx3G@*(9f#{h8T$;W}PHkU!Dq$dbvYb5m$-=(myn$r$o*GClkdq7g zk_(R)Q<4Rz;s=tm<>c%qr}QOH9Y~%oCr=+po+~HM?Mt4AB%{;i(K%fyXD7*Ni;v%V zEp>Fe=hU?0Yp%IdPaZva^!1Fa;FP63n=h>F%?p<8l1ruF;;(jFUg zW=v=Hnem->$z$jArOgexR4;_KOnpu$OH)CzrfuG`J-?M?w$0^Pn_hgjLsu{xWd0XI@NuV1Jeo1Ed3ERT~ zm}Z7cK%Fv*U}msNNg!Mt4nhxYkRHsOC#My5d!Ma+rnV>lk_R3YSR}usOkMz|3E}^* z?pwp!y0Swh3n3wd(EH6(0`W4KH~4|g+h86x7=#~zU=WxX13rLp(llkJ^FvIg4LIqC zAk%BvonFh+^b79S-jM0E!PDLeY3J*l&0#EALt3UoXYT!Ozwg?PliW;y^saS|B*e>M zrk(a)KC#X@`~6sZ?X}n5do9}TASwqfb^K?y;kEs=;T9Zh-`<29Yh#4O0EeKQ`z1^v zf+J|!REGCFS|h=+V&9d}a7!%Idl#TGyIb?cYE5j0`ez&;iBRM)62s|0DF+QGs6;b^ zsE*RofN~ND9r_E^co3ZSM~IPU*^MJZEeLyEk%#@eO=YLRGjYdKZa|*ga{Uav!AYCD z=25k|%E*_xGTN^HS z;hML&j=(R!ZuyM2@dZ54)7b?&*Xd3NY3w6;eMHqqOZ#??tWeYbn~;g6z8_-q>^^Zt zhmdI})s=_vBuH|1*Cw_-*e`a)CpD-_Is=N5P#5PBqAH`MWxuXUv*hwQmEIgdJILhb z2nx{1FBw5xC_G+7z4vRR2`u4zPK`FZaVNo3g!(5$?Naw-K3K9442Nf!sHI)2B?6RU zc`y&W?3dK-mj~M;9KSVUE3dXkIKlRCArdPLoVMQNNxlqKKs^S_>`(y%Zz1uz|D;zZ8lNnpqko4 z>j9rZccn++xkhUb5oInMa98BPp(`LrTrAXPRZeRU5@j|Vc~|7Xk$1Hfj=U@2czfh- zvae9*o(!zF_h{%XN80xcBiL7`?YDFJt#F<wd#`GAkR;woTh4jftu`GsSD zv`}#fWVj)}a2yAzZVnv&2>FG_ixojUIZthQq7W-Aj|GzF;(y1#ow)H^x8YX^A7FvW zaxj~y@EiiqGWRUg*Lr#<5n`EpmUM4YpuHEETLY$Hcb9P89b91E2v8bswoZD_a^JFW zc}akp`{|+vlGjL7jr18fRDowVch2G=`Z4)yXumzi-46N@YG)w&+uB(V#d;ox%@X}k zV|VSGQdy$E*^I%Fwl~q=5jGoJ{Q}qO!9jM?)r%LrX*EcpvVw!`vRi;RhiPpVQ69uW zmK_{qm)#LzgdK=}4iJ4H>@cK&unV7y*#8b4vqN&SN5@2LC&a;yxfh&`Hzh1=KRss@ zY?8s<*FT$UJGQ|Ci=Z!A+l2)!-KAF=SS6p^NZ&xvmxu?7_rft{k1)Go-EL{^z`Ee>HRXrq2yUiGaoP+v z`Q95UccHMG8!9N+3v6k!taY&SxY>4Z)_P00PQ+$8Q*Cf1-O8TJMi$PaBgR}t zoj1*2Q|V?8-#qN5{s8STUA9kEru~aJ;i; zrpJ{Y5Dyx2i5BN&P5yeoCoRAy8V?aI@`)Pc6ZN38UZ61maw?R+URu6e=k3G;xN{MX z2wQ}*L6X5z@CodeFC!V73_MFRcG7Zfh;?EKL2&J}ha6z{z=ONejbDm||Nn-X@5SCm zhO@VuUGrsgtt!OX)_0XRRe@k@XBXif6t~8DsSC{jOyK2imbcygrt4n8CwA?oT{zq4 zk@wvuT^o{j+a`J6oir!@zqgtaZXaoTOya`rTfcYn`sc_n?Voi_d{<_@vR!7rT+t!0uO&g&+y306dptJM?;(Uc8@v6K*J?y-5;Ya zVe)MV@B|dqfnXjJu;be~0B|FKz{$9O!|>M*UH6 zbPJU8!yzuXj14&qz)rIO+}MLmuxB?s)C1*qG!FN|y*WWEc(AOf7?Y129~!te?(j#B z4Ws`5O${oXl1?{<5uWa*y07}>O~XS?j($*xaoy44=pG#94nl33HTf|FW=%IxcoUM! z!MI{Lxk`b}fKjRv-4f7E0EaTE3jD-&;%?_w;ad%YS_G_gThQJOCeXrits5G1^gF0g z>Z=gOFR8vULV+R<(9gk%#QzXoBLl4@hRXD@4ILd&t^uasfw$bE7sf9Me4?R=Z(zp% z6#-g-qW&1){u?IAf7u^_9pI0042@rPP|)ap*``_L252{bP=iN6f0jQsEYq%@;UT}E zySJCqIPwq^UPFMEdHqW06=+fFWss~5LgHWrgfajQY}d%ueFpQ$?6~F@jbEXAz7c7zf~wz7BI8o~?!!z>lHON6huE z?!m#(rxXJQ3+9Kfq0KKD7c}GxQ6~a7&oNo}918o4H$}WhsH;P z=Q8UlFDCV4Oe%t<0U>3felSQBGEL);tQZ*Tre0<(iuj|$Y^(UCu#&@L7)ac@{X)QD z7^O*AF~tshAk@ELZak$& z&p`-w73UX}54}u%12cUS0iGWPLNhqXo@3=Q4n#o`#4ZC_VJLoWo{~e(fD0Uyl#Y&U z^}f!PFz5}72kz)BCsUx_#~KIk28C@j=3N={M`J^<;UH1=ye<22ohwanEy7hem_BkT z{w+HWA<->Yu0VPaYRl#w;A(zYN!d^@`Gq6Yz|a_lW~RAJgDX_7C)i9w<=8xevB`t` z^oo=>QxA=va>L@h#@Zk7WAbZ^- z>~U^?bRh7hm!ZccFzIx04u0WvHi1Z>nGn-2x*ot?yim`Y{J|nNvG?mkpLX_Ae&`d# z^?q3}I7@|_wtd_E@-6-t1Ml^-NpJ(CX}g2nuh<#_Bk{Yq>;4C>H~9}-xBCT9Prm>r zehN0tVP3rCm$OtGZkZ1C4zR80mveUN!(4#6BN;?eARZay-l!9h09B5l5kV7zDFA-y zOCchQ3|WBSvYZX_9aI!tN~R5xt+{ZZqi4rpHB)2KTb^B#gf zLGb4QfJp@I5X?aU&W`#cu<)CDn*Fh1o@|eafdzujf)4&E67Xm)D>k-20{#34CPskQ z9HCY)2Os7jK{a9=&q4h?x;oJ{IG+c25D0QxL>fROj=~9Gt;#7; z`2<74=T*N9nkG0BfEVz`;y2gHtmReE^A1C1jEctTK#o9%ARcRRpcneJ*AD;Xeo;55 zCL8J+VqakH<6gjJBIv@spmD08OE*=DX`Mh&ji3#|3kdoVppr3a62aF2_=P~SpihMZ z>?0Bl|KTfQ)&Sn{|rHBLWdE>neWBuOX&4gY(2ar;Fn-W zB8`~nrM{0|e*?fD#r7V0f!)n6w!0!hD;#!PBn)>fm0*<_30#q_j{fN{1j9K@oyd0RBiMc2I$G#N`TfZbHz3 z;6((N5exzF3%jp(4-A5Zy8!))fs>l3m7QbRy$fMJ|q`SdPW zO8qXT_7us>zeBKq;0}Vj2>u4aI|zP`;KvBQkKj)c z`~`v)1P>8>55do|BaR~|LGT9%($M2y5Ie5o8+-Gx58wI`pn^Z@FVWYA;A;rd(IW#v zCW0&k2N7f=$U$%j0efuow=qT|etv*&&k&&fBkCD|(Z}$K!rjfk7dKz`^9FZa{*(Vs z{IY-CS!y8L1A{=Fa3xIrIUqRRI~aAs(W8H5|KRcsNQ%XJS1uzmugUbZ;V2_R!(+Y} z6|J_Cn3SgnKvUimfkYPhxtb@6{z~s}_dV9|gwfCRJe^^7=;n~S-BWwJi_sQN);Np5 z((Kd6yC&vhz4{EI&+vl$TAv3Vr6fY*lc}feQ^Vgna94Uq3d-ja3y7>>veqY2ypb`L z;guv3Ng^Xjo~-m~^r$jswri&ATdj;H-BU?4ISZEe((b3dmvcYouPsbo4bjvd-$#DW?q6HLmLE<4lyrT??Z0sDqzJ@+68W-RlNW`YVim zCgtgjZ#pKM;X2Zr<8*|2-Iz!x6%k`G9ij6ntcymHQnegKQqB-Xi&tSM3Ol1XH(9?S zjrB_7R-|!^G-0y(A7v?%wV>JYjsCCqyLFyQPvPwZBFmku{#xyZ6wao*QfJD%sx+bk zjVKv#9h0bX8C5U?2clF)c+qH{=rL#m)55Dzl&-c|DW(RH#xP^B`w?!V!ySjV2g_QSf_ial(s7qhl z0ZBBk?I+N>-sCP=h**StD< zRq2D`$y2n_zOK;G`b?t8f=E*B^eb*KLVj?;Ky-(R^vGnjFGlB*xDOM3_ChAnA0shk zlQlk>-eq*x!y$e++-4{vvg4Dr8|q|tA)RrAQ6HUbT;o0tGwLJop|!f(=&YlR_Sj_e zhQ{Kqp|g%LnzG5JHSXgmqd5j2x)gUGD9d7W#{q#hAd;QVKF(-Q09J7yWsLSXeCSi% z6ZD}IjJ}+X(5x%9lZG zG9y8wO|Akkgb<3C`A*Bb?f2RlMZvz}t!op#+Ek)V^=i|mYMhnR3NXU#O1Wj5v%Pr{ zPC&K2YrkjzBk5x9lJG;(1JR!!eI(bfD>bvW8Jo-UP3qS5x$BHF)mga~t9WDL>l4$X zZ@fD7Dumg9JvGal_0DQvg4J35M(tGXboul(SJHGXk($@#CRh244DN1ZWz$kK2iNt6 zTf#Y^tN(48v(Bf}&yLTGyV~FUP3K9tXUY0pVe}cBUgsHIn)BpFplM1uv-3WAqE~JM zCWo4&DV8X z7Bp7#RI{rbP$J0%jb-(1cUzTOpVjesy7VUop33^Hm(Pn&x^-^;+@Gea@f$e?U4=mpLZUR5qC8Qir zBlyDOpuq;N&kTKi$klesK4*73?)KmA_l&h7bs}rWNb7Pf zi;g9X+&Wq36UkRZ=5=#|FFqN6D(#afpa7O#I34m@yoAywf0Z3Qw9QDb~ z0fQbYcgO8GMxF!65i6g18Bob5)p14?+&8QA8kkJbWdeOy8}FG`jfLMf-B#n3 zxVLlfnO2O2jH<{fUNe|c39siQF%&w5-;n*pmcwWd-MImcYPos+_6=t=oYi}jp3CU- zzO8qcyDI0k_w>#fTA%kQHJj1o+_~&lxKua$ZeMmvXid(KwMkq-*W!}h`LwN=SW1|< zQs*fbf2wIM!3rC`31FO;NjU6mbR_|bubGqJ*uoW(oX42+o%Pe#UT*}=54vP1wbtM= zdJP%Gkii(T7|lVa#HUa80=d&?da4+Gt~17`PI7m`ezG%aYrcP3N^Hf;ImB|FiMxP# zzuxpjfQd%(;u3Gw&DTLc4c#7EGSThlS9|)rJre|e<{M7os@SjwSTZm_;FBueXqswr z|K?)v5AwaIy2z<6P!pFR!7wIB$eeJibiUNxcDM6(=aPgz+p*f+>+K#V@H1V7`x zadzV7L1#O;X%;d`R{ zqD9RIy7zTpm)iEgMxVY&^e;e!7>rPf5megQ?3>x)6YxUed&T#Qm$XEA5<-gBwR*sD zg?3i=mJZMuRXo3i#zk~Vv#-s(=Dr3dB>EWL!RZK}O2?=yK6Qduol4ZHo(NC(g2Mgk zLv;aMAP0LO0-n+Udq+`fuUJpSdZ_E8+@n7-9{Us1?S9TcD$Zv*8fqA=CDF+~O(Kkp3NVfW z_3V7B6Y43~c;zWXp2Df22lD3CRbiKZr~%l`a?cb;;o}|M_HFk3&|cSF%_;Ejs(6Td7LpV=1{<2-OQd<5<>pk5+zSJ=r0Bb>yf^`_a^Zoan5oH0yvFk7qw zD;}48%RUeMcJ930#rL${6VE2i*(d91MefrGp-}Qn4#;@xwfWaPFL84eXH|zjbD1Oz zU<&39APoKU{kvqsPE4CbD83KFn}GUx;}5OGA7bM3SO;H>@(uep?7)xh9H^+kJ%k#W zM4jn{W8mr>Kn<@xhv;(__zSsU?nkj8A^IbL99~l{G35r0Y{F$7)3GI?_rz&(;xyi7 zJKKi3WY2Vv6X)oQFOU;25Yr2QDTYq|Gr(mn>>O*1UO02Ok9;zT35*MXfk-qjU@BO+ zMocA(h2A4)$Pv`ji{{IK&$w_AZ@M1Y%;z)A2XGA#4u(t|0Hpy2k4*E}$kUi+n`WBa zMT|PbQw)=iR}HY3QI~qvM~V8VS6xBW6<&2cQP(r-24J1JWz4sh&ZHz2?-Zp*4W4Yb&#KRVI{rn zwh$Q2yrtx2{sY_FNub`YwZq$bp0u9-Xz*4NNiTE4BJ%>&(Kx-;R2&9-@P8GnOK83B8@{8F|J;)FXh|r#z-xoH#Jc9XKJ3- z=KhU^>cu0BrWTf`TAN3@;8=t;Mcrg`h@!2@RX53B zoSQv>i!GZk^Vm4fj4nS*Ol^p@`rz~eqpxYc$y2;g{a)?;TDqd0nC%$A5Dbuvs3^zmkr%o3l34UT-QB+R&XwC;RMVGcbRcV}Yxl|HToj#F0_%x-K3D*$388ln z>I-rW5NVnx#hY=2WE@$nTv9BRFV>Na8lYVDd^@m>ziR9zFc|KWC%Bb##$m|23G)WY zLttn*kd<>?A<_&_wKwY+$vU>^fO{1D#Vbr!{U@5>*m@ivLg=*^`H$pA zR^XWrEVnduzlD(tKGS-lSG{29{1W{UOT*4xRkV5!q}>w z=%86Xjps-WjK=Nj$3db$#g6q)1I+C6!h9A81Iy&lG#e1N%80R<9k{`CU~q?~=!L7q zP{k62=lcfJ50FGoVjN3MCo!)6vxrDwx0SptxxYq?HEh1m_x(H&H$Y2SW{NJYBW9NM zF@Ox;VD=f~xgSV?Vj?H2VYL>Y;WJr%Mo@1Y!MaF!u-|Zs#3WAEz?yHn^=9TiQuHJ9nD4@U*+c|M%=7rOqXA)z8H;O3?Mg zYVvWx44DqQFrWh{GaJc&GO&a=@wf|M;_+o2LO0HcCn(Q?&P<1J3zo%*`xzwj*wzSu z8&4D{T~*^kbYgyDFZs}FOPjZ)kH9aXABPMIP40^Bh<1uuxh3psYrD6#pTIAXWv}Db z5|V$Z^brHbSqzv!30ziehoq*$HR_Nj_X-$`@sIbeaZ1*|AT(| z>}6tRm)3^O9LD(|6@u~&%@1)K0US57IR3%J`xDE;ABsK_(e3>to~;SHKG|Q$byhy? zgO%s`ne%jFF)Ri^2pX*+(Mg^RZ^mKx!BoDMoZ^1z?(m)A1+Y3(LYFp?)Mh66G_cm| z$rc)5Wuv#Ul~lH_R<@zln##8K3l?)1yWhIR4vI63*5=jbu4r=^ZT`ZU6>SlzY(wU= zeZy+ITYIN=r%-XVt;5?kK;UP62?btB8_6l&Ru9&gD01(y$u%$SmAh~lQeXZgIy2RFLGibT&1J< zFA=)fmT|ZLPXEr~y?VOUd%BmvFO?O9TGk-s#GQ$qV){*Khz&%tt77ZZh{!0}Q#k;T zw1o4rTcs`iB>6H|Z3t6#r;MAH77ITpe!qCB`or1>wRFn>v9N8G2;@G_<~3&#bJjz1 zHZ0YuW>Tkfr@MWi8WNO3y3R7O$!-S|n~8!*y&3k7UY&n+K|vDoi9R2;bz=2y2@{)v zzHv5B6p725YI4?0*KCz5+E#*E!dSW8!4>sE5SnTiwXiQo)F)WcB5W3Bcye*aQ0P5W zK@L?wgZC}>(UoV(p*EsJ`J$#Bgm40d62&3bC$YS^9FT;>Wf4W4Gjcu7a_i#U#sAVd z)d>4x`X8Gt(;{Dd$~y<%m3~`FA8sURO|T@Nj`peI-JtLM5K0x3GCikAQpt>LI?@%1 zqQ&@c#)D{4r*{kAeC?4wku7ZI$`%Pv?us5|M(vC3bY3OVSFwU5yj}*@&_J40?JcMw z1vN|8mQ(4L4!X936r3Zbb0|1EkA(+;j)-SW1F4SL3aApF^~?-g3$y2Lp7Yoi3w|K* zp1erlr@1l+Bo#q4N;x~v1wo%#T5MNk_GP|9!1CB;4WfB!EeLyH9SIx?9%6#IO zH$+n+ZsU^MCnRaTz^^G_a{pT0r(#}=Vp{h)$T(1nQ+u=4RUIyo|8AxC=xG8!HOO1O z279t!9tP*eHa^HMf=pESNT93Rh@PFCG#j`pdDG^yZO!k4`tQ^eB}zEhS++4Or7hxq zy4Ftgu(J&-w10vsq~H#xSDCt^OobVSfjPYfWhhNy6w%I6!@n5+@2$()zq$angdcP< z$J^+O6p0_jgiKiUvZs-J;jvv?c(rtpsBpKwVGa0%ZK9ji(kny-l9sPEY~OEz?Yr*^ zyd@0;emYi|q(-4v_?FDCEeOw;f2xPG9gsa{MKNuPgL$lp$6-c@n6`WP?u;51*R&$tIsE738775VJ4^ zp*(T2&aW(*JcYtY(Y(lWWU-s*Pq3)Dv28PLhfLmE+p$u+WzUZ7l)!dqxX^Sfx8SXA zU#V_i?YiXc8YS>EjA1-2^tUPAYs^|PW--QWbV7-yH!f!-E{BPOV=y4eWaUo36O^R? zUU19n_i)P#=fdBIO_sIfgL9Ft1Fj06ISzNWVNHI322nPUN(PH`HwWTCjmLH&x zoZ`eqKRVLxGbee?IV^8Kc8@fpx93B8HsZSYAS`tCR`)qtn={$$SEWxj zJklrAAkVwX=&Rt2&?8Vwnv}Ot$7stYo7c7R)8{9f{K~Y+Q;&2>?t=?)+=-({ngn-@ z=VeAyI@$C&bVB6;sKH(gIo;Sp;#nD^>GKG&Fz%TaPX%n_`NJfgnBWOaa39@@37R3n z0RO-W5jX-8CR<&XD&MG|s;717j5GsIX^EsOB0a4B-D3+?OX}rb`a&0dd5|ctIz@iH z%^8El49*&0uj2Wjwpsg(o!j*JGb^4UT&iLj*N&Al9oWzkWJUGR3*6kX&uB)7_6-?y zkIp8(HOm6Jp_82KqBSo%YkWE*Z32DohZZb!Q9a2!wcJX>xrfF+THEie{ij$ZQC7TP z#3(Dgu@xk?0%Vk4xrN9r->SIVbf;;d5X1{!YvqcylCf5=%4xKpmTICdRq9g z&NTbVo3FT2xLtPNp~C5EuiCnzwsO1LZaYyQ^cfP}syXv?WH5B97y3da>`JX8J?6HjM*njwWXLws;|cBjZ;m?XgQ3OGFQG=dAc7jl8c$++_v(zq^> zN*H1JS+6n>BG8$O>JW)J+mDv3>9>rMvfh&_fb?G2J*!_ZSpmCp8GjwQ`7w?OXlE zUeb*psJwMA68I%_vqgm|CaNE+1#8%$J44~yC&80DS4iq0hP7b!$06uW>c&m@^r;KP zd=Y~h_FvLqqI=;h-S7g5V^Q1qDWIebBtu)ah_;*=a*!l;p$$D74MEM%snsw^9SPN{ zYu`wKOLyh|O27r(N4Ku446}!44%6}65xS4XG0w~1i1kIQR-*AFpQm6kV(Bc1`(QJR z9vmUsYtCA~In&v=0cjiV7~tFw&cs_EADa?xmCTjEfTK~E zGl?RU~7O`ZbAe2AbE94mG;6MdGa zXCa4VSAgQ>SW?wN>&{I_ti{_rF%RR5W@@G@T?&|dfF!0*S7I>^zkS%aGwgi{MG0Mo1u4BkxfDyRUuwwOuAUI;$CG#6#?im<4m2!Gg?}#%D^y zr+z(JKH2ca(>fbZ)xrPblN$$krlK#N0(X=9g&s6rs(;IxAGVhD0 zhk3C%Uwm?r&of?x2v(jJd+9ElhI>k){M67k@KIFfyy%!luchzE5%_=cxl8dEuv zBY9wtlX1@+aEI@U5^ygRKe%WDcYGcMcX$pagZmM9=7T#tuL{Bai1>R=>G_q{_>um| z{&8yf(!c;-3hJMq!@qmSUx|g=_@HBUG-qI7G<#rhq(=yK3I2)T1uk6hKlxLDphB15 z%0CSH4!UbV(SI<^t~#8OL1R)OTwB3E)@(oo_Z|*gF!gPJVCx4I@w9@L%C4C&A?Wb0G2xr7? z;uf?)z3q?fa|Nc=k+y69aD7*@s~_8UknS$++NuOh;sd6MTa4SXgwGb#D6;J@Xy$3p zIH)ya+uwzxpu?t}tYD2H9g&s&PvTuj4u%u)x}g<@qAsKaU0PLbQD98tz-GJzUMqws zkHYxeTOYrVcTAAL>*J@ZoBTq&`Z)SOQFF*yFv=}qZ=(*vJ$bn2eRTko7dZyI`r*=| zXzViO={ciDvU0&B6R%kpw;UzUO&=JIq+t!*l`2&c61G2@hbxPf;o>-Z~|Yj67>-zQ}c?& zUf$AF*Hn8ru?8+WI(iewh7%_!P-&7_SJ#@@2Tq3*(^UaIE3Ir+HDKpsn1Z|5Pyx_P zG71`2M%W9ULbfRWNUlVF4czO6i>;tdWq1e`I`j_U8x;vn>z51+^*Opn2hiMyX7jz7 zY8ztj-GaftNUk1!;ngc}kJvBeRKS3(5l|oEj|oZVJm-^ScQ*$mreL*1b3WC@#7gnf(sDhA5lH3FvfD1z<+sH0yLp!$Bv zW&`>KpjBZM$;jqC!N3hXHu1m>yq`iMb1+GaqB0A(dMD$FRnz^9C=PB;c@@bkiez^& zqsVkd`oyv~WK%M)I1Ub0t;eVN5>nxaqc6dN``xl<0#TF%Bmr_QUW)-y_?YumnzgIG zA&K!yOe+$T>l7nN!)1?h`Z6Pm_gN0U5j_`NXI=RFn*8u*c)@C4S#$6`s;ai3I;_ffo3xYUs3u`a4jv(|VsJdW=hnalAR@J~ z4sdMF>XqgYY0j!NH_Rmi4)45^@44|u$5y4K!0t#@lTFUz$)^AFOvIBLKHhPW%qvCZ=IhbQ$)wQrS&1XP*;SxCqy`iXIbHzON;q%*LK^N$rC~f)Arl)p z3OrjK)PQ3&5S;oaKqo4DODO7R;09_CF*(kEGnL|9LP}W00&FbmB>Gn)s6)W2WU3Fn zVHsXgje!;>$A!MZ;qI}lME3ClRP#K9dIbxQL=MK})e*M7Y!XmdRv*L$7WOM!lKX{X z7HZp9*tA;Nk=(+uW$Ajqq^qlM z95h%ty1FO`&f~O*08JdbJUDQPjqMkLjw4W5#cFB-Dl^pEJ<4gXihS8Lt*YB1-})KuL{DKPVpmZ(j& zL#zrD_GT`H#Ni%iii0B^;-f&$#s?i`6l}!tkRxPoGCvmZ`TUQ0$v@>qf6UYWl&AbL zPxUd+^f52tV_wF(DS_tcpA|*$mCt0ceBCo^1OS_qW{{oWNJ$nEE literal 0 HcmV?d00001 diff --git a/src/optimization/__pycache__/run_optimization.cpython-313.pyc b/src/optimization/__pycache__/run_optimization.cpython-313.pyc new file mode 100644 index 0000000000000000000000000000000000000000..0b42911d0ba8d529eb1065cfeff37656b2355036 GIT binary patch literal 34471 zcmdtL3v?UTl^|FkKE#I~kN^l01VMn{|Njw1Qj{oxU-2i3EEYkDk|2--KN1O21yCQ_ zPBWcOhl+dFP~1L(JWd29o)*-}jNx=<%*^gi=ydNy{n_m*ff1lGla|x&?e*!idn8LI zj{EGH-utRh07!wd>`A82*%DRt>V5CKzjxnz_Z6O`rpht6{`S9pd;ai64EuZZ!e0_D zAb4_IfMIW7IED-Qu|5GQFrsi$e-fF*K82(Zp2B`npO_R!(j;Uun&08F%XOFCc6Dd9-Lvw2&g4Ligq zP8}(FwgistthX_V%ft3S(y!`Mlj=SVsp-oivmj1^C+o0tB2tT|;8He}ji=%=Hk5Np zhszI)0|ZilAg2K1du7CClV+33(mqC2mkAh!-RWF%+uid^PS-Jb7ag~|=in*6XkS=y ze2nl^CS(+oEAtD}w#B7s#{#2RSh7#s-1f-@hl^3TmKH8LxRj{~ye!mexC=2%@Z=s# zzljLOF`(lF#B&#cjfnVg?8)Ob*8jNDNqj0`KP5u zQoWEWf>f6OcBVt>*CAC5sVO^C&>N(ZkSc-Hw1lZPNKM|+L*=vi#`len5y|^lZsPmn zloC>0bl{IMY#Q)5Sd3wDy%rO~Gz^yzBDg0LVz^5Q3EX8wGTh}*Qjv}|METl?5lMJT zRVu)eG;TWZzeMUhzN5}mzGWj8`!|az;BC;2anEl=o^M8;|DJo!Vr6N^R+U4O=Nwb+ z`HK!C%Pj|13yv*0t5zN4l5v?_nqGmlC8u%PF*WZ({ADaXUv_{*Ofvg#{_p?j`}DsC z+E$Z{##NEAs?E46G?rT7wG;m1tHQDa5K~<{Ln7*~>Ws$f>T2WCihFs*ZFG|>&MCXw zF^wY2WQ=&3oOilOB$|wHY4RK+=w$?bjKInWIvD}Z2+pHS9P>&S?U!xj(j^y@%$9OY zGlF?W0A(&Q0tmYp;RW{cSY9IC zHqzl*S#Z0^Jg5UnfoltHYuL6(M#=SWIecBWdvEmy4b@ws8a56E!fEMiJ>k^M*Q|eJ zT^EEivaW0YSi5csYjoGGe{A&#!a0VU+PAbGOE|CKru8kWSFmj`Lrp!mdi<7vX+z*2 zqzexO4Xs{bIM?7)fBV3;-sEfdUGPo$)qcAl4+sO6Kto`1BV$9jQNLl|z&C}P9e0JB zEmXnami|~+Z}wfdG2rb9=NJDxzkb_TL>0Gh89RPqEPj%X8H%0=F`Wr+S$S~F%=tu! zrD>i(zp%tYf2IM$&R}OUCjXIy98XI9Q!Jag#{j1Y{VLc6Y!s6vP(-~#xU)RJVx zCJ;b|&Fp$3;lpEndX`5Eny8;1uHDz$I)*R0Bdhx`%q`5gY}oEV>{) zp@M+)pFnyhq$gAnkp8ETu7dQ0N&?c~fb^A77oNc>obx1%+oiaSdDUYjt>^Zfkev#UpianuB`$F}|Ch~bojQ_X_coR?Oso|94 z>$~u(M0KANIa3Jnh>0&3dvdfUl19w2I6lPj6_4}L6+%KtE|C#%R!iP-J+4&`+arAud#d;hJHAGZjy+5jfM}Eeqb#;i-kE zasv2(yGCH(+o31kG`O{ho4GRIUd z&{?h|FdWuE2KLk_kiG!vstF{<>_x=z2O40ahL3*%xlW)bh)g1)M|e!|&o5MS0u8ly z&3I0Al1(`1)a08y-k`xEeUb&4sa`SaJ#Z?*!;wHSHR)CvZ z1kdGUHC{ldM_|<1eYFZ;QlJ4C)kHetXjdk9#mjhp zb!rK1qMnT$ezp|UUOp_~OA^|J=m-OKSf0%v#~JOP|lC&DZAw+5slY7N9n zMKPJ6CZ~?jCBjiW4xKu%C3E-E@kf94JWKKD-y%ElMaEdCp3oy}dJuX*=+XF#7gh_M zc^o7#KD+@nx&GBGO2co}qBOj~D0VDHK|YkwI?X5Y2?Mf6h`do3(pZDnIjuwFCbD{n zJf3G^qDx1{rjNL1oct%YBXaq&Bf~sQu_vd?qB#=rx^8dD`&=$G0;fLDdNIX8pGbS* zF-Z9kVM?S)LEjfd^*z`%&H|zUuf^;32|=R+o{>g(B{>WDo)V@wy_)@87|mY-M)M29 zSpOUt&2caqp9^E*mw>VGg<))X4vd9yFy=fL#-c9{BfbxBtcfn45k;}KcFgYf^H!9z z7_?pqkqVl}!SX;Bgp3d>C6o|Cc0d_Xl2F15h*F|#UvxYKIg8i7H}vE%g~GsB~E&$SweFJm=4v&Mi{SHyVR z7g}McA*vI}bY!hSPt@RttE01NkZj-|Bol+@BF|SM5>_3DD=@&~Xq|Xut>&rftc}Z? znEz^b&BEawjr&CHh%H8EJ_PhW60=jCwU8@;H5;?2aSLnl+WGRaCvM`#ti8xgy5URUqahJwUxGu9J_jESF?`SDR64M-J-3i7? zkO`RG^70FoMZ+X${ypk#NW@DoSZ`uT6Uyt0)fh*k#5Tpgh$T@4<544PuWiGKRe&s< zXY94cgn0aXdu@+$&)92w<=btq?Uiqry$16RR&I!@L|li2zz#MUu~%(u**KaP;4IgD z$=cYmJ8EOK+}^dZD`9t^10P#HR^#9=BV%XdZ^=4sBpPA#xpQ6C26-5ZIUkq6t6N zl9Yrwn?WAu`IOj`XvT-|<7g)mat^S&c`$Mk%_Cs(?7r|JWDPi52qDzmLa3v2sRM)v zV${R|e0bNtk!mT9J-dc{ zu*GsYJ{q$94}=n|Iv&bV1m|)LK0Y+z;zKw=w4!xWo~j6I3kUerK`4?2ePjt~hAG^k;d!(%EC}9{G9< z7CuDuz}`2=EV>eoxRrcu04?+zaeJ&SsPh|K+4wuZip`e)zjmB=6K} zhkhGc&3X)W$|vws_-WW}ISz7sCMw6VT{bzQ0+$oG8*+lIrf`^X4nqnoq=8nj!5YSh zb)W#V4#E3b+}*AqiR-4 z4voAXFQHu6+KC%jKG$$hdx)afAQS#eG8j z{1;grvyP9t(YntdOg#{70D=>VSI_y?FZ&S;vmzJ)lL*GW)(fLhf=^c5Ivx+@!WV+_ zL=jk}ykOYSv5xy9~j&wKek;CrpanpH*IE|MfM_W6r`Nw>kaHsVQ+KjiJ zMk~HCzp5btc!k;E|Jl!>G|bLN-PRNAtQ7VxttURWH}-_RS^7e~nI!Cc^k$OnjXgnc z&RS0HBLc!53IH=o5gP}7b7t=n z_{CTW=j>kf;+Nub%~_{LzPlSk7{=@A9>Fw>50EFpy6ecld$t^YIf~hN!0d~bx)LpQ z4pw^4AM%oVopLJn@2t^Fm(ooKzlI&GZ-e*X&Wwd`2Cfmt>timq({o?vWa zXK1j315kc}6`UTy&>MTqb|P~4xm>*ubZ;5{doYiJ`I8a$Eb+Aj@(G-@=X}G-r{G^@ z=Pa`XH-q_o_HB;mb| zY*gG3$F+a(B3#RHV;t{-mT_LN55E9& zbDyg%z|LfOjnfuiMv@3(FRXOp>1bwSo!^C*Mk!Z~^Xm*55SY>o%nigQ#4D=%(YD9hFYAU;WY_u#1PBGa@htut1 zq>jt;F1KyzJfpB(g5!BhwrMybH3em(7gA&@bR-W=HQ&G{q6s6QMR#AYT8&m1-~^6S zKnNzmB96^mi0Eit2q$t1u`fyo2UG^D=jXg=m6gV!&pQWlkxUe(@+1sjPdCH-Ts|apzv??v* zwfS$(dxgFgI<0D57|t`^u>H}B2Y+>yk~Rj+a3c+nb1Jez2tf>!!G9JTKLQ}&WKz-; zoUE9L;tNhHCBVw_jw?Eo;z~z%8sV6sh_^B9s0Gk_+yw_AW8ownMu?uBtKdQ;Z#+lC zOhf>?q2dVyJ9`0?pt*NnyygVW`}dqSSe~5#8-crsZFd91~kIS3P!{?O6FfakuSEMc(pN@-p_~2BI zU<7Q)-51Eaq&|s|aPti~WfR5d`1eO()ZKkefGxpYd>e_*T#ZCw;d95%h>S#Acs4B# zT5%#j)e}{~N-oYA>jE_LT z4Eh2%8y#AOmX|VE?_z*I#;0)yRiiS^w|fzLLkw*L&mwDl7HbmD2<`}uA^BZ!I2obK zO+q-Sx|R_*cL=iV;H<+50!+3shxU?V`-*$XHNWbxU0Nc~yXG7YcTC2k#p|IDE~Fbi z1|g?RBR{+nFnzD$ZS&hxe$)46Zddp!et6|GmlQJmzNAF(N4X_}kCDv6QPLXY$0+*a zyKpO0GNJ|0jf{BSIqkU2h@JLDIJqo@lc1oooGXhC5;P@~Jhim2vgmX%ayYH(bh+)$ zDF-8VuPiS(7}2?FIDa~LC<2Ytjy>@N3}Wd}T$mmSV&sBs?90!LBHQkWz- zVVX2O1C^mRZBVz1k)gAvm*%J4bBt_$`m&9KLk?+NESyP2VMtpem^I(J`RcXeRBXZt|sCjK`?p##V11IInY^^@05zN-ov8-J+sTi#arhW{k% zT@`I?+|umd8aqLqJRKZ6v);X(mGhR_JM(WVd>4XQW$Qg*jrO|rn$;537HSk=0idA;gGq5Hg^QgUBRqws30XPl#)lM2ahJ~&=6g0HmPdWTpf#uadO2L#mN}BW8Cko88--X2Lry5O{MxJ;KIS5Dt zNfRLZLxhf&3nyV{iqrfs{KTYg#}RlD&YqeeglM&fmpxvS5Rh5|ZO4th3D-DrAzvb@ z$#G~uVcwDG=m8Z+E<83~i6>yJ=Pp8wi*PY(!NX{~01-aB>=G` z-5(=rlu1aiIx(^vZcOr=-DP*XNwN(^jiDQI$SaeQM^Unb6*~A>!Ge7f>U4oA;b0V; zWO6T3`}DaL7aS>PWbUP9u5A!ykVG+ChE+2Bd2mkFvbiD=(Y4Jikv1+XQ@f|cv5h=i z0%6nzvsol}s+=1pPs_SAtjMAj1%GmRTL&Wew%Ir5zi{h(pr0;294t5z)LGV#dWNs| zZ!0oA_16zvJ3!^vg%tI)qCT9hyD7gR4`nyd*$n|pFnj-c@3x`DFWbo2GPGWkdxT+C z_Rm!%+l3{5!4Dezm+!o?@e+M_n64QPR*VD-anI4$2Oo>E{DY4qm@<1kowb}8!5AZ- zTydUvkTy7l&q(`h?ZmPEu8|#s4`M8K@EM0dnwXAv5hc3*^GR5+SE znlvij`Avk#5M+uUg2);u4ZI+@d`C*65rRZ^0f`2C5_X9(5*a0=4&9K6NWwf1i;!=P zy;$;*FfGOy(_WN{4naexG;H=@zhTFq;)#3s$-DMjtAV2TO5QEmFl`KP3O4ID?Hg5e z{n1cKFICbDIX-(Lf+AllDHr^!;fHhp1mrmcJ{k#)#@(GK&qMSuN=HK)*S*{m_lC#s zT5k0DRzf9BbV<{OV54rszL~LU*+`{JxP5LpDbf2AIha|nRX>mN8L z$CJ-Ie($*TZr!au|H?a;?pz97;0YLbg|6ue6?RjF-H?Op-M?1HoG;QbBw;ROkZ^=6 znz}bJNsU==oeO9}br!nLvT52JzALy}ch|mIMb{k<6%SLz!;mRLx&2UdM+Xs;faD0C zG((_G@&-0425A*0lvWHRRWd3V744MPLX4G7mNI3@jAWj*a6U&ekvSdDn~W4@iI|^QK|Y)y1Xlx-%VwAb6PC<*wWOT8_X#A zh+}?g?vj0B8q7DVV_|ja5*Tf=^Q3*oeaS(dXT+Bvr-Ko??X!#+X3{H*OtNcg4hp#$ zrPDq&H$&QI9Wc#yI-xq{!n|X~0UrgKuu z`Wgo=msBPoVUY-9GGhIRzpdp_(VoD*5E;~e`#RQaAh>Z$H zP+->;o47$v)*#|*fp_*}6#J6k(;ty3d7>JSslum&K_NLtu0tE-6?9ugw?9NTRKXn) zle9by^MEC{)6NLN*kG-QJ4uXyWKu5LrX7oRCw%S5zRXQ)*0=#*dT_zMI5}-^+nX`@ z9AAZ2TxZ}0jQLS|#%q_qdD$zVbd_{=RVclhPOp9>#8T6!G$Wm2dW0pVjtU-The(UqxqEhtg~4^qL6zO>{~D zLO=1s=)XdZ;K7VhgnoiUKa&}ibcd4{<*xt$AZ5rJ1))WP8P|====Ei|!4!IB5#$^R zBb%3zu&ZA%546wD&cH{G*xo`7NEErABfVXg#O0RH%r8VNWT+)kE3cs30`g7v0j4Z6 zsToO+116Df24iH$ctytK4DcEZOzc_65LH@tgyb}oXOz**HcnxZ--c*DRTWE(=vW8I zW?&?ovMUplKFCjg7jDoaR3-Tqx}jS5i0|x5`- z6%B%Wo2-)2qdO?ti7gzd8P|Zlp^Nn z>3hT}>i95gk5;~0xp9FT_BU4P`o2&}KULBX89sBd)@X5=0EobHUreSD*D#UsXwrdJ zJOvk{!|HT3pu12+6a@k5EOTSmg7~HKiM|emJvL&4ea&gFi`k3 zynk$g8$Xj?(FF5?(Tv> zcV}07XDtM}n!B638&Qm8B<*YsjHr9b>1Gu05kKg+ zvm5jSBf{sw$PvSL`sSw@NdyJ~lRPuqff_H9FcPj!M#41-nb;QLoox`kqxRsPYfeD^ zuc(hA_uRY_@CTFfD+>+)LbTwRfmJqgein=*AYhSuN=^$r~RP%> znjTdqSADt3hhCbgr>Ih6S zcWJC*z?5qZGO-%#5tdw&_NyEte8>)c(eB&mi+0Zj)Edu4?*-p5r7H0^zSDH4iB{D= zQD7O(*W}*Dn@u;Ge09DVzdc}~%#C#J{*bbXQZ_x7!WZfu%OKZd2?`}c=))W%EM;Ug zGdaf3*d^1iS)`U$!TafzaNofZa)6NKi2nluJ2>Jb%Ms!|Iij282=SgA0ceB|QNmJ0Sz-Z&%K*2rabteTAi)f5^I_N!`cC*sCi;_}ELkpN7p}R%23y~5OA1|U3-QJ`$_#}RT z-%e?316}X+yxT)-TM;oW0+yPWZd~$RMAY0Erb^oA{KFyj5lVd|N>Ga+LQ#t#M$%?W ziQp%taQk=~VAw&|=z`*(K*|ofBB=u{x+h(c)Fr{9qG&vmmAVA_I!d*e)FsHp*O7o` zyQD6jR!N*XjTf{SnSKGCA$O4{B4`Pya~x4;JVqU{G3!^RQm-X!D!#4#cUlPO3w+&v z3zb_JAl^Im?kPIg@+29{=@HzJ`x+0_8CjSAY{hx6A_+Q{6PoUQi(d}1o1B0xWd;C}Q z5Rv~oO8Faf`)hPVnjUN#wo0&pxEiex`7t5}T1g=t?iJEWmawnKK4wS28uTI}NqBG@ zg8zgeP8gFv`fUz?#Yed?%Q^!pK~(x5HJX zV2iUmtS)35qL2DAT%?^O?N?mvP7b5!h|K3;+Xjqv^oiX_c{UuGGK*bv^D}PpKe4Fk z;vG>K1qv=$Bl7^j5b{xhIWTMUyVzYKs1rcqOj9O?RqAWX6sdRy1W? z6uZM63BuO27=+P%h2=_+w5 z%ZM+dSki@BkaLW9-nIC5nO zW^k9X@ug5^__#A8gD1D_$7Db-x(4W4UYK_?DQ@^g`hv^0ykG}3rcT>|pXXTGlE7&1~E3a z_ppKlzzdPWMLj_~bE%gcj`OHdHUoJUwUZJ-F3a@jWA5xe9t%C6#Dc}hB1z!5z$dWE zfysy@0RAr8xTqA4RFO9ctcN473S|_w2(+9;rm$jXm7d?wXtixgwo+k=|0Axucv z%r4T7Uc{5EuLWpNK)(T5LWz;=YgxwOa#!E+@x z05ZvJtZfpOwO1D5Awh+Kj*?3*8w438df$OR1&WEFjgi74nT^#o@>$XWE17IK88ik$ zHR#%kARF#i9`Z(xgFlpYdg(W0!7 zNKcFOzR^G@g+eXc8QETqFVzpr_4R9gPen=M{ZFL|vF@o(CeD7O#&Yv0k(T|)`@i8o z^8N<@k@q+FkG!)}Vr{m!-=+otAcy7zT=Eq5>6 z?W8TkboFqsaO8&UqiE*YJF@{h${v`7?7_mrUfFh0rJo241dA-8qWMkx-!xJ9ICXNG znwbms&(q!Wo@3kQ;!yDfT|5CPA^arlVF!yRg60X&z^^5kzWk{W%h88&Ptdt1sNu2D z@Hp%N2XjwsWuJJMSM~^|YWhcVtfNpWPZiW5Spm{uK9QTU1 z3yMNTQ*_Z3HE}jHF$ueo!J?^P!IU?7yR^#hxU;lzAy{hhrhb|v(DVqNV(9)*T@q4P z)9UJnW#!?b;;^|eTvT)W%&jwl;SE8cJ6O~jE-4MywuEaBZnqxU)YFy|)WoS^>*?@5 z%Y%}dfGk*Y;6ZtP;9#))$b;f)|K(tD^W&V<9NnkFbWPUN5?CsL&SmAnEi>ma`mFsY z#{^hr*;-#XL+uf-^)l+FwLvB=Z>{H})b!U*ee)Ehs|u!8hf*8p)P`Vc;~MxUWnXuG z%jp~Q_uR4tvmj=@>8nGJ(y&bJb?vvbPZh~($tQBCe%-=;*&h-8sREOyt`EPO<;hrA zY-i=Z)p4`;&0e1a`AP*cVCzRs7w-@1n}S)*+gklw_{|eQUWJujQ(!ny2z&Fv zyyl>`<*`_pCI18y%H)qx6X3igQ^QgSe4-yUyEgDB6Dugat-Gb8D%<`d>1VQ^%BbUG zq2s4$_&ahsXgagjO)2tdks%D(%5RlZRY(4!;Adq&Eu%)pLnCKs_-nTX3(lf!23nLK zHkI7g-_k!#63OI1^R)CwsaQcZtm_m=_?jLmF>~4N%3GB?GMQ^2lbOwg>H)aq%GdI) z=26Cd)c&JXe(!o7CF={xhG^MPP&T~Q9!}2jWg^If$rWLf**8d=>O-dev}yl_g);4@ zOl|8)2z068!-AS%R>ya%yu*RIjijH*e=L9Z@bwF~s@8_6tPXHvkttqFyPC#iR)jT1 zUm2yY466#l>YPV1Os0K=iDmLnQvjH&SsuZwIgipYu>ye0(qBJw?Mz5hPiyLb*b>wn zTI>FUW8qZQYiF*W*;Zy=Pra5JQkK)oa(^FP*}geM!9sk+>7epVSfzd>#?thUFku?> z797c-dHGfOJylCkdVp#jq@~CHPL_66^Q!jmWXh{rPi|0_$A-K*I;UpCvaY3MZCp~; zb^SHH--Kc+A1N?#8UQ9!tk$gF~Y zjBjN}RWeMR1^lVjU2neD98#HSmDz8iRAx$58B*2Ks@kBcKBU@5tM;u8gjM+=RSB&s z3CQW%ZVDfTRXIvUtPMO=XRQq;_{|Dc-u1R?ZLmaV_Z0_KW$T^cbk)_DA1JfKdZWKC zFidAQt#^bq`CFO__+bl`+ADm$c{|UTfWtS6ep3G9a^%ObIZ7Ysqf7dNdHvVAKZ<6) zbnBA;0ul)SB~Bnb-P`(nPLe{!`$2>_7mvUtt z0wPj%p%mDQHTqQdQ;MGC@S5PM5z#K?n<<{tL8y;t<&GEx=|IMza`Q>R2PKQCA7LEs4iO@+}0az4&E3H=^JT%W4Ng^ z)HHCfX&@*agnPiWZlY34{RPlJx@mxx4n7ek zQhCiA5PRspv<+yRmJO;%oOWkvD`WUimwc|j(QcYO9bT1p3K062S~|1f&eFY%VOli2 z)*BY3eu7CdC66(QShCjpREec2Upw*56J807=obSoZ4Oe}5$eQAIu!`>kkbRugsq8Yr*QTt9T}P)Jott13TANqg7ey;i{m75X0L?RNKP|vgwAPBGYqLpa+M%7*>?s(aO=!E9o^H{5IqHTTlZy>}~u z&0}F@@nd09OWHrCsjfOamRILMBd2CjS^MaeM&PWLb^+CXmTtC(n&;@|Icok~u=)IY zFPHuh@nsRMEMjyG;k;7sU|3&Il-Ov+i7^fS?fupt>CW6XS z_ZlYrc8|$p^G$-t&FhjQ+tOHbT>+`@e09Y_&8e3>%%Vs&< zK2DvrQ_dyI*g_eXsg%6!((-pw@1%xGEp(|R+&U0y9i>}GL#+!>937>Ew7&d3})m#&BRK}Z@+Tul~BnMy5vZ>sUy_X zPdD|4noiM8r>N6swwi1~_1Sw(wt&fF_vV9|(tw(>(dx4ROlhVkWlLEISkmTrM%TN; zDVeVgULB-#6`ZyVw7+-s-J^FiC{Ucb{;d=MHd_D8DmR)%!dw65AuzH|A`CaKbvE#-l*(in~G-^ip&4{j+B#YIw;J$E5?U`sg|R%t{GK1OR`_2I)pXD`9YIw`T-HGPCgkkgQg=O@ z^XG|j-dA@$&cU*bkMvk--j9d2jJE$g6u`HvXQ<=2&+K#gJHD4sThH9e|MAd0qYd;H z;3u4-2Tf~A15FD)P$_yCB9--lQz00-8heB3z0qz`r7RKB%lxN8Rp6g*38q`(65F?w z9pUFCKFlrQWRhw=3MUg-akQ=vJQ9h;1s_Or9%{|&DNm)CJmodRHw|8ef8d_1Wn(NX zmA|IAs_=CCF0U(SX>C}c^&Si=O2V1h*Q>5o`J{Aa*?PwVO^ILV9|>w|*SjBRs$hTv zW_nmq6wWUQ8%*fW1QV4~^jB5;I760^2@}9P3OEuB}pQFynPvQN@=Q^CAx zkDST3ds04hTmZ#j%48fDs8e>&QO*hqs5cv_QGzZVyRWv=jtk%jmSF%FFjFQlg2Cj| zhfEE$sUc))p-nBpyaOJ2*jO1d?xT(SLdJu%@t`O5DeNxidtbSqQ{}O2=N0%4+|R3q z-Og;&&GZ}Tus2rcpA3}W&pre)Zm0oBb=5FU$}0=$YiWILNN=I_;F}xNcOqLf*CW|3 z09@TK*nh1LEE5gjDWoo<)kQuBFqgVEq;8_sO&h}d>Vr^~xd|Gt+8>yCZ{gjAKXcx% z>;kNKx^5iucKP;&btQiBkA#82U~M~H+D_{_JRJzWL-+IQyn+u*MX*y;SQ{#ArVE=x zg}roP?_E={aNzDRT{z_J3hRq*4&4|61g~yb?tEoqWpnDjzVC?`D=0=CsfSQix8Tz> zEMq_rRQ5moof5pO+tV9T*U;)3=ui7b#=G5rrl8da@2lHc%717m^OtNi(}qrug2|uo zq&zq}8Z;1pm9klP<0XI7hGi?i?JrY23d%q{&>MZzH-->P=0C{E^L2yi|3GJeDS)o> zLAKtvpUy6SkelzDp>wMq%fu}ErM(~{gki}4Z4osQN9C8^{K}25!0g8UZV?Fk9TNsQxwb7I}#LL*CRr5kvS1mRt!#hS8T5OeqBq0j=qg4$I4jb6V&iOry6n zq%hM8^9Q;@j|5KAysh^5PpJ5sxm1h z=(nJ<4C>4+0^rqIa016r8ZuPVhU$=^_Uf_qZqIHA!v&S0f_l23K2)&pTAD}VMQ!>t zx3#ylK~*U<7tXE>W!KW#wSlY+RWSRIM+`^VKE+^vJqCad0iezF2mw0>EpV3Pe{OlJ z%V+irZ5u-20O(M(I)YC0eKL9f{F98JYBFL(JvSgfEN7ginmwbmmwSPOl&mK zRqdOWt&(nP;tZuO3@fv*r(a8dSX2=%uG}sx^P7Jpd{6d6*+%5#9yM&GPS~hn8-3tx zu+~nO+o`GPV4-8XsNCQ1Bh`C3Kg`)^47HBXts@i;>v47px6`eY!MZ8BVv2Ii1dC?3 z&82?T4_xnD`O_;KnxEwVIR9=I)Oy!Iw~PmCPS9oSDxUfD3uX%!Rfj7NvO~$ILPbXA z(<;o|2*k{Z|S{r!Hkl1Nm!cpTKd&=uY{^-zgx7PPD_uoL;MyD@tRyT zIj77AXiY0Hq&DyROV?gnKg#6ygProA_e4;4@?P(WO^dfVlv_>bRs(Zc{J(cU_b}ai zV!fZzoeXR9L)vm$3#&pKGH4*E?N~n=&M}2@D(Rd`|GAAas;V=X)3x3QtaAO}^@E;+ zVV&`2^^NL~u9nu-26BVCL!M-`$Y!`@2$`B_QxnzP7c}*I(!x3er05O63Wma)OE;G8 z8)~j~!wBgLYjxgMT3h083Toj17}yXu(!4SlXMN{^JH5$aU6F78N9Mr!V0}AX92*+y z&_P=N*gL20oZ2wmu?34^tY|}Hnc-&2jTGOoKjRhwBZ3Q=y&Hl%FKu+(&prYw*l-vu zZGFD?B|5jlKNTpyvlPrdw5i&hfc0o<)Jly{(WBEr%QR%oF+I#H;TMPjy#v<Ovs zTN;0IKoyuk2UVK3)L^da;*!!fgw(-fbo0=C_3_8KSbpteJ(gAKmjoo6upHM7XFZ{h z#b6F(`e3B3;%%s`kkUkhAVdHCg8RyH7`SqfnE*?9@l@pz{|jZ)ugv`*QxCewf_m${ z?y=2`8xwx*hG{GBAl*IYkx_c<1JIb4Hl*NF;gLSbD}~fjI(u1y$UM($o=dt=viT}0M>ddMg&N!ns#x} zrhwL260*@QUfRt;CFTBBbl%Inaqh=8!IIv)I5m8V9z0E*oua0f=(Eeg!Q~sx zk1*_#z>);B|8`J^*ZZ&a`z(IJt*$_cr=LL17@eF1W{9rj2<%2HZc`U91J zBw)HT_WfL1*?3>sirA+5)0e<9KS0S^K6^q)VR7TLPmTyMZKvR~Pc90u+6j0{u*Nxf zNHN(l!Dmm(vCI*8JC=lHoQ9_k%Qy!QJ*KfkVm)R$EBMXRax6vj1QUo0!dla-DUXB@ z`hPz;o&?t4XP=xEVEJQEVnK-IqT1WAym7&&1OnpcvX*~#q0RTd61OK?97&SDEmd0_ zQW1qo(BlJRy2X(y-L8j6ru3I0cxa^`s&gz3z4V`&;H^OV|6uSakv@{bqf+{7Go;l? ze{%?mH({S;cja4V48s4`1P{67U$T0{==J^X+ziW1y5#-t%3&=VF)qxr%%~;rj~lxi z*a({na^y(fw-t72*a$~cie)BQ`2LJUz&+(^Hh#t+x6D)~PeB)=7xD)KXL!Mr<5FCP%h^wm zC9YpFf#78?__g!Y=hl`Mmp2_(J~wYhLb#(pyosxxTUtGVd=fdRG z?3yI_$-LYaF^d)C=dqrskf`|BaZ{on3&!>EITzss z+I?lO@w$<0J{)^OD`;1C*9EcqT~Y8F`;=DfgBDP$F+on?Gjt)h)Uv!|tEQ0^r_ntJ zF5g#{7MG0izoutI0<3|MYcRl_Q#S zEG(}|tE$-FvIlqk*%fkWa(=ZMN55@P8qwA}RD#+ux|bsGmPSe+Fmk!es=ebY$D$b*)3Z0HvRkKKhm>Mp&nvHE?U8eXSTSrJ^tBv@yfGe4C@`@3>S)A?~^iIl1dr54^xEu`v3=~QB^J1mliL|Pg=5Txqu)QoMVn$i>nmBr!IX0%h4 z{TLIAB``SJT4U~QKDt9)V@86XR#Sn;nem$GH7`)DuT+& zaO%-LvZ=OH0TOjMwQi3LdbS}o+OMDXC{TN_RSFQpMh+N9DVZUKo>u6C3d36GN2%Ff z_07B+dBM~Y*pf-ndL^N{ZW{hlde*u=5~X{@$mI8R-X6L&6c~LE+-bK&?H`F?av9Q< z)4KAXiz>LI#8MkL_RKqNHuw+)(z^eg0pGc*?f@r7#HoXDK4O7XBMxAB1v!Z_qGL-= zMbC!+CL@TylbE4i_Mi9}DMm3g9Ejw%L$pXGB{d7^V0~hiznN zXq0tsWkjs!p9H>#;~?FPE?!Ajvq~DgyyoE$b%fiI1hCvr=AjUBoniTsf##4Fwb|}q z+@V!8lrYkRP{+!Gqm8@;IpGg;eG_hAml$zcq&g9lsqlGBoiDFJawxY1eK2#1ehpgeKjb|c_5O0y&F~r JMLHJN{}0?MlSKdk literal 0 HcmV?d00001 diff --git a/src/optimization/model_builder.py b/src/optimization/model_builder.py new file mode 100644 index 0000000..4d3434b --- /dev/null +++ b/src/optimization/model_builder.py @@ -0,0 +1,1479 @@ +from __future__ import annotations + +import math +from pathlib import Path + +import numpy as np +import pandas as pd +import pyomo.environ as pyo +from pyomo.environ import value +from pyomo.opt import TerminationCondition + +DEFAULT_STEP_SIZE_TONNES = 1000 + +# ----------------------------------------------------------------------------- +# OVERVIEW +# ----------------------------------------------------------------------------- +# This module contains the full Pyomo optimization model that was previously +# implemented in the notebook. The workflow is: +# 1) load_tables(...) loads all parquet inputs from a directory. +# 2) build_model(...) constructs a ConcreteModel with: +# - Sets: sources (I), consumers (J), weeks (W), days (D), shifts (S), months (M) +# - Parameters: demands, tolerances, capacities, mix targets +# - Decision variable: k (integer steps of 1,000 t), with x as derived flow +# - Constraints (in order of creation): +# 1) Daily delivery tolerances (power plants) +# 2) Weekly/monthly plant tolerances + total system tolerances +# 3) Veredlung day/week/month tolerances (Nochten/Welzow/Total) +# 4) Mix bounds (hard min/max) and target band (soft) +# 5) Deviation split & max-deviation pattern controls +# 6) Veredlung deviation splits (Nochten/Welzow) +# 7) Smoothness constraints across shifts +# 8) Route exclusions and capacity limits (monthly, loading, availability) +# 9) Transport caps from zugdurchlass (if table provided) +# - Objective: penalizes deviations & smoothness, applies route incentives +# 3) solve_model(...) runs the solver with configured limits. +# ----------------------------------------------------------------------------- + + +def load_tables(data_dir: Path) -> dict[str, pd.DataFrame]: + tables: dict[str, pd.DataFrame] = {} + for parquet_file in sorted(data_dir.glob("*.parquet")): + tables[parquet_file.stem] = pd.read_parquet(parquet_file) + return tables + + +def _match_kraftwerk_row(df: pd.DataFrame, j_key: str) -> pd.Series: + if j_key == "J": + mask = df["kraftwerk"].str.contains("JW") + elif j_key == "SP": + mask = df["kraftwerk"].str.contains("SP") + elif j_key == "B3": + mask = df["kraftwerk"].str.contains("BW3") + elif j_key == "B4": + mask = df["kraftwerk"].str.contains("BW4") + else: + mask = np.zeros(len(df), dtype=bool) + + row = df.loc[mask] + if row.empty: + raise ValueError(f"Keine Bounds Zeile für Kraftwerk {j_key} gefunden") + return row.iloc[0] + + +def _match_ver_row(df: pd.DataFrame, label_exact: str) -> pd.Series: + row = df[df["kohlesorte"] == label_exact] + if row.empty: + raise ValueError(f"Keine Veredlung-Bounds für {label_exact}") + return row.iloc[0] + + +def _match_bounds_row(df: pd.DataFrame, j_key: str, zeitraum: str) -> pd.Series: + label_map = { + "J": "Jänschwalde", + "SP": "Schwarze Pumpe", + "B3": "Boxberg Werk 3", + "B4": "Boxberg Werk 4", + "ALL": "Gesamt", + } + label = label_map[j_key] + row = df[(df["zeitraum"] == zeitraum) & df["kraftwerk"].str.contains(label)] + if row.empty: + raise ValueError(f"Keine Bounds-Zeile für {j_key} ({zeitraum})") + return row.iloc[0] + + +def _tol_from_row( + row: pd.Series, + d_nom: float, + step: int = DEFAULT_STEP_SIZE_TONNES, + min_steps: int = 1, +) -> tuple[float, float]: + lowers: list[float] = [] + uppers: list[float] = [] + if not pd.isna(row.get("minus")): + lowers.append(float(row["minus"])) + if not pd.isna(row.get("plus")): + uppers.append(float(row["plus"])) + if not pd.isna(row.get("minus_pct")): + lowers.append(float(row["minus_pct"]) * d_nom) + if not pd.isna(row.get("plus_pct")): + uppers.append(float(row["plus_pct"]) * d_nom) + if not lowers or not uppers: + raise ValueError("Toleranz nicht definiert") + + a_min = max(lowers) + a_max = min(uppers) + a_min = min(a_min, -min_steps * step) + a_max = max(a_max, min_steps * step) + return a_min, a_max + + +def _mining_week_from_date(datum: pd.Timestamp) -> int: + shifted = datum + pd.Timedelta(days=2) + iso = shifted.isocalendar() + return int(iso.year) * 100 + int(iso.week) + + +def build_model( + tables: dict[str, pd.DataFrame], step_size_tonnes: int = DEFAULT_STEP_SIZE_TONNES +) -> pyo.ConcreteModel: + if step_size_tonnes <= 0: + raise ValueError("step_size_tonnes must be positive") + bedarf = tables["rohkohlebedarf"][ + [ + "datum", + "JW", + "SP", + "BW3", + "BW4", + "Veredel_Nochtener", + "Veredel_Welzower", + ] + ].copy() + + bedarf["datum"] = pd.to_datetime(bedarf["datum"]) + bedarf["weekday"] = bedarf["datum"].dt.weekday + shifted = bedarf["datum"] + pd.Timedelta(days=2) + iso = shifted.dt.isocalendar() + bedarf["week"] = iso.year.astype(int) * 100 + iso.week.astype(int) + + weekday_map = {5: "Sa", 6: "So", 0: "Mo", 1: "Di", 2: "Mi", 3: "Do", 4: "Fr"} + bedarf["day"] = bedarf["weekday"].map(weekday_map) + + wd_to_date = ( + bedarf.sort_values("datum") + .drop_duplicates(subset=["week", "day"]) + .set_index(["week", "day"])["datum"] + .to_dict() + ) + ih_dates_welzow = set() + ih_dates_boxberg = set() + + bounds_power_plants = tables["bounds_power_plants"] + bounds_day = bounds_power_plants[bounds_power_plants["zeitraum"] == "pro Tag"].copy() + + model = pyo.ConcreteModel() + + J_POWER = ["J", "SP", "B3", "B4"] + model.J = pyo.Set(initialize=J_POWER + ["V"]) + + weeks = sorted(bedarf["week"].unique().tolist()) + model.W = pyo.Set(initialize=weeks) + + model.D = pyo.Set(initialize=["Sa", "So", "Mo", "Di", "Mi", "Do", "Fr"]) + model.S = pyo.Set(initialize=["F", "S", "N"]) + model.I = pyo.Set(initialize=["Reichwalde", "Nochten", "Welzow"]) + + I_W = ["Welzow"] + I_N = ["Nochten"] + + model.d = pyo.Param(model.J, model.W, model.D, initialize=0.0, mutable=True, within=pyo.NonNegativeReals) + model.dV_N = pyo.Param(model.W, model.D, initialize=0.0, mutable=True) + model.dV_W = pyo.Param(model.W, model.D, initialize=0.0, mutable=True) + model.a_min_day = pyo.Param(model.J, model.W, model.D, initialize=0.0, mutable=True) + model.a_max_day = pyo.Param(model.J, model.W, model.D, initialize=0.0, mutable=True) + + column_map = {"JW": "J", "SP": "SP", "BW3": "B3", "BW4": "B4"} + bounds_by_j = {j: _match_kraftwerk_row(bounds_day, j) for j in model.J if j != "V"} + + ver_bounds = tables["veredelung_bounds"] + v_day = ver_bounds[ver_bounds["zeitraum"] == "pro Tag"] + v_week = ver_bounds[ver_bounds["zeitraum"] == "pro Woche"] + v_month = ver_bounds[ver_bounds["zeitraum"] == "pro Monat"] + + ver_row_day_N = _match_ver_row(v_day, "Nochtener-Kohle") + ver_row_day_W = _match_ver_row(v_day, "Welzower-Kohle") + ver_row_day_ALL = _match_ver_row(v_day, "Gesamt (NO+WZ)") + + ver_row_week_N = _match_ver_row(v_week, "Nochtener-Kohle") + ver_row_week_W = _match_ver_row(v_week, "Welzower-Kohle") + ver_row_week_ALL = _match_ver_row(v_week, "Gesamt (NO+WZ)") + + ver_row_month_N = _match_ver_row(v_month, "Nochtener-Kohle") + ver_row_month_W = _match_ver_row(v_month, "Welzower-Kohle") + ver_row_month_ALL = _match_ver_row(v_month, "Gesamt (NO+WZ)") + + for _, row in bedarf.iterrows(): + w = int(row["week"]) + d = row["day"] + + if (w in model.W) and (d in model.D): + vN = row.get("Veredel_Nochtener", np.nan) + if not pd.isna(vN): + model.dV_N[w, d] = float(vN) + vW = row.get("Veredel_Welzower", np.nan) + if not pd.isna(vW): + model.dV_W[w, d] = float(vW) + + for col, j in column_map.items(): + if not ((j in model.J) and (w in model.W) and (d in model.D)): + continue + + val = row.get(col, np.nan) + if pd.isna(val): + continue + + d_nom = float(val) + model.d[j, w, d] = d_nom + + b = bounds_by_j[j] + lower_candidates: list[float] = [] + upper_candidates: list[float] = [] + + minus_abs = b.get("minus") + plus_abs = b.get("plus") + minus_pct = b.get("minus_pct") + plus_pct = b.get("plus_pct") + + if not pd.isna(minus_abs): + lower_candidates.append(float(minus_abs)) + if not pd.isna(plus_abs): + upper_candidates.append(float(plus_abs)) + if not pd.isna(minus_pct): + lower_candidates.append(float(minus_pct) * d_nom) + if not pd.isna(plus_pct): + upper_candidates.append(float(plus_pct) * d_nom) + + if not lower_candidates or not upper_candidates: + raise ValueError(f"Keine gültigen Toleranzen für {j}, Woche {w}, Tag {d}") + + a_min = max(lower_candidates) + a_max = min(upper_candidates) + + model.a_min_day[j, w, d] = a_min + model.a_max_day[j, w, d] = a_max + + model.step_size_tonnes = pyo.Param(initialize=float(step_size_tonnes), mutable=False) + model.k = pyo.Var(model.I, model.J, model.W, model.D, model.S, domain=pyo.NonNegativeIntegers) + + def x_rule(model, i, j, w, d, s): + return model.step_size_tonnes * model.k[i, j, w, d, s] + + model.x = pyo.Expression(model.I, model.J, model.W, model.D, model.S, rule=x_rule) + + bunker_df = tables.get("bunker") + if bunker_df is not None: + j_bunker_map = { + "Jänschwalde": "J", + "Schwarze Pumpe": "SP", + "Boxberg Werk 3": "B3", + "ISP": "V", + } + bunker_rows = [] + for _, row in bunker_df.iterrows(): + j = j_bunker_map.get(row["anlage"]) + if j is None: + continue + bunker_rows.append((j, row)) + + J_BUNKER = sorted({j for j, _ in bunker_rows}) + model.J_BUNKER = pyo.Set(initialize=J_BUNKER) + + bunker_init = {j: 0.0 for j in J_BUNKER} + bunker_target = {j: 0.0 for j in J_BUNKER} + bunker_max = {j: 0.0 for j in J_BUNKER} + for j, row in bunker_rows: + bunker_init[j] = float(row["anfang_mo_di_t"]) # placeholder, refined below + bunker_target[j] = float(row["zielbestand_t"]) + bunker_max[j] = float(row["maximal_t"]) + + # Determine initial date for the month to pick correct anfang_* column. + first_date = min(wd_to_date.values()) + first_weekday = first_date.strftime("%a") + use_mo_di = first_weekday in {"Mon", "Tue"} + for j, row in bunker_rows: + bunker_init[j] = float(row["anfang_mo_di_t" if use_mo_di else "anfang_rest_t"]) + + model.bunker_init = pyo.Param(model.J_BUNKER, initialize=bunker_init, mutable=True, within=pyo.NonNegativeReals) + model.bunker_target = pyo.Param( + model.J_BUNKER, initialize=bunker_target, mutable=True, within=pyo.NonNegativeReals + ) + model.bunker_max = pyo.Param(model.J_BUNKER, initialize=bunker_max, mutable=True, within=pyo.NonNegativeReals) + + model.bunker = pyo.Var(model.I, model.J_BUNKER, model.W, model.D, domain=pyo.NonNegativeReals) + model.bunker_out = pyo.Var(model.I, model.J_BUNKER, model.W, model.D, model.S, domain=pyo.NonNegativeReals) + + def bunker_total_rule(m, j, w, d): + return sum(m.bunker[i, j, w, d] for i in m.I) + + model.bunker_total = pyo.Expression(model.J_BUNKER, model.W, model.D, rule=bunker_total_rule) + + def out_rule(m, i, j, w, d, s): + if j in m.J_BUNKER: + return m.bunker_out[i, j, w, d, s] + return m.x[i, j, w, d, s] + + model.out = pyo.Expression(model.I, model.J, model.W, model.D, model.S, rule=out_rule) + + def y_shift_rule(model, j, w, d, s): + return sum(model.out[i, j, w, d, s] for i in model.I) + + model.y_shift = pyo.Expression(model.J, model.W, model.D, model.S, rule=y_shift_rule) + + def y_rule(model, j, w, d): + return sum(model.y_shift[j, w, d, s] for s in model.S) + + model.y = pyo.Expression(model.J, model.W, model.D, rule=y_rule) + + time_points = [] + for w in model.W: + for d in model.D: + date = wd_to_date.get((w, d)) + if date is not None: + time_points.append((date, w, d)) + time_points.sort() + time_points_day = [(w, d) for _, w, d in time_points] + prev_map = {} + for idx, (w, d) in enumerate(time_points_day): + if idx == 0: + continue + prev_map[(w, d)] = time_points_day[idx - 1] + + model.bunker_balance = pyo.ConstraintList() + model.bunker_init_con = pyo.ConstraintList() + + first_point = time_points_day[0] if time_points_day else None + if first_point is not None: + w0, d0 = first_point + for j in model.J_BUNKER: + model.bunker_init_con.add( + sum(model.bunker[i, j, w0, d0] for i in model.I) + == model.bunker_init[j] + + sum(model.x[i, j, w0, d0, s] for i in model.I for s in model.S) + - sum(model.bunker_out[i, j, w0, d0, s] for i in model.I for s in model.S) + ) + + for j in model.J_BUNKER: + for w, d in time_points_day: + prev = prev_map.get((w, d)) + if prev is None: + continue + wp, dp = prev + for i in model.I: + model.bunker_balance.add( + model.bunker[i, j, w, d] + == model.bunker[i, j, wp, dp] + + sum(model.x[i, j, w, d, s] for s in model.S) + - sum(model.bunker_out[i, j, w, d, s] for s in model.S) + ) + else: + def out_rule(m, i, j, w, d, s): + return m.x[i, j, w, d, s] + + model.out = pyo.Expression(model.I, model.J, model.W, model.D, model.S, rule=out_rule) + + def y_shift_rule(model, j, w, d, s): + return sum(model.x[i, j, w, d, s] for i in model.I) + + model.y_shift = pyo.Expression(model.J, model.W, model.D, model.S, rule=y_shift_rule) + + def y_rule(model, j, w, d): + return sum(model.x[i, j, w, d, s] for i in model.I for s in model.S) + + model.y = pyo.Expression(model.J, model.W, model.D, rule=y_rule) + + def y_delivery_shift_rule(model, j, w, d, s): + return sum(model.x[i, j, w, d, s] for i in model.I) + + model.y_delivery_shift = pyo.Expression(model.J, model.W, model.D, model.S, rule=y_delivery_shift_rule) + + def y_delivery_rule(model, j, w, d): + return sum(model.y_delivery_shift[j, w, d, s] for s in model.S) + + model.y_delivery = pyo.Expression(model.J, model.W, model.D, rule=y_delivery_rule) + + J_delivery = [j for j in model.J if j != "V"] + + def delivery_tolerance_rule(model, j, w, d): + return ( + model.d[j, w, d] + model.a_min_day[j, w, d], + model.y_delivery[j, w, d], + model.d[j, w, d] + model.a_max_day[j, w, d], + ) + + # Constraint: Daily delivery must stay within demand tolerance for power plants. + model.delivery_tolerance = pyo.Constraint(J_delivery, model.W, model.D, rule=delivery_tolerance_rule) + + model.shift_balance_dev = pyo.Var(model.J, model.W, model.D, model.S, domain=pyo.NonNegativeReals) + model.shift_balance_abs = pyo.ConstraintList() + for j in model.J: + for w in model.W: + for d in model.D: + for s in model.S: + model.shift_balance_abs.add( + model.shift_balance_dev[j, w, d, s] + >= model.y_delivery_shift[j, w, d, s] - model.y_delivery[j, w, d] / 3 + ) + model.shift_balance_abs.add( + model.shift_balance_dev[j, w, d, s] + >= model.y_delivery[j, w, d] / 3 - model.y_delivery_shift[j, w, d, s] + ) + + week_to_month = bedarf.groupby("week")["datum"].min().dt.month.to_dict() + model.M = pyo.Set(initialize=sorted(set(week_to_month.values()))) + J_power = [j for j in model.J if j != "V"] + week_tolerance_debug = 0 # extra weekly slack for debugging infeasibility + + model.d_week = pyo.Param(model.J, model.W, initialize=0.0, mutable=True) + model.d_month = pyo.Param(model.J, model.M, initialize=0.0, mutable=True) + + model.a_min_week = pyo.Param(model.J, model.W, initialize=0.0, mutable=True) + model.a_max_week = pyo.Param(model.J, model.W, initialize=0.0, mutable=True) + model.a_min_month = pyo.Param(model.J, model.M, initialize=0.0, mutable=True) + model.a_max_month = pyo.Param(model.J, model.M, initialize=0.0, mutable=True) + + model.d_week_total = pyo.Param(model.W, initialize=0.0, mutable=True) + model.d_month_total = pyo.Param(model.M, initialize=0.0, mutable=True) + model.a_min_week_total = pyo.Param(model.W, initialize=0.0, mutable=True) + model.a_max_week_total = pyo.Param(model.W, initialize=0.0, mutable=True) + model.a_min_month_total = pyo.Param(model.M, initialize=0.0, mutable=True) + model.a_max_month_total = pyo.Param(model.M, initialize=0.0, mutable=True) + + def y_week_rule(model, j, w): + return sum(model.y_delivery[j, w, d] for d in model.D) + + model.y_week = pyo.Expression(model.J, model.W, rule=y_week_rule) + + def y_month_rule(model, j, m): + weeks_in_m = [w for w, mm in week_to_month.items() if mm == m and w in model.W] + return sum(model.y_week[j, w] for w in weeks_in_m) + + model.y_month = pyo.Expression(model.J, model.M, rule=y_month_rule) + + model.y_week_total = pyo.Expression(model.W, rule=lambda m, w: sum(m.y_week[j, w] for j in J_power)) + model.y_month_total = pyo.Expression(model.M, rule=lambda m, mm: sum(m.y_month[j, mm] for j in J_power)) + + bounds_week = bounds_power_plants[bounds_power_plants["zeitraum"] == "pro Woche"] + bounds_month = bounds_power_plants[bounds_power_plants["zeitraum"] == "pro Monat"] + + for j in J_power: + for w in model.W: + d_w = sum(model.d[j, w, d] for d in model.D) + model.d_week[j, w] = d_w + row = _match_bounds_row(bounds_week, j, "pro Woche") + a_min, a_max = _tol_from_row(row, d_w, step=step_size_tonnes) + model.a_min_week[j, w] = a_min - week_tolerance_debug + model.a_max_week[j, w] = a_max + week_tolerance_debug + + for j in J_power: + for m in model.M: + weeks_in_m = [w for w, mm in week_to_month.items() if mm == m and w in model.W] + d_m = sum(model.d_week[j, w] for w in weeks_in_m) + model.d_month[j, m] = d_m + row = _match_bounds_row(bounds_month, j, "pro Monat") + a_min, a_max = _tol_from_row(row, d_m, step=step_size_tonnes) + model.a_min_month[j, m] = a_min + model.a_max_month[j, m] = a_max + + row_week_total = _match_bounds_row(bounds_week, "ALL", "pro Woche") + row_month_total = _match_bounds_row(bounds_month, "ALL", "pro Monat") + for w in model.W: + d_tot = sum(model.d_week[j, w] for j in J_power) + model.d_week_total[w] = d_tot + a_min, a_max = _tol_from_row(row_week_total, d_tot, step=step_size_tonnes) + model.a_min_week_total[w] = a_min - week_tolerance_debug + model.a_max_week_total[w] = a_max + week_tolerance_debug + + for m in model.M: + weeks_in_m = [w for w, mm in week_to_month.items() if mm == m and w in model.W] + d_tot = sum(model.d_month[j, m] for j in J_power) + model.d_month_total[m] = d_tot + a_min, a_max = _tol_from_row(row_month_total, d_tot, step=step_size_tonnes) + model.a_min_month_total[m] = a_min + model.a_max_month_total[m] = a_max + + def week_tolerance_rule(model, j, w): + return ( + model.d_week[j, w] + model.a_min_week[j, w], + model.y_week[j, w], + model.d_week[j, w] + model.a_max_week[j, w], + ) + + # Constraint: Weekly deliveries per plant stay within weekly tolerance band. + model.week_tolerance = pyo.Constraint(J_power, model.W, rule=week_tolerance_rule) + + def month_tolerance_rule(model, j, m): + return ( + model.d_month[j, m] + model.a_min_month[j, m], + model.y_month[j, m], + model.d_month[j, m] + model.a_max_month[j, m], + ) + + # Constraint: Monthly deliveries per plant stay within monthly tolerance band. + model.month_tolerance = pyo.Constraint(J_power, model.M, rule=month_tolerance_rule) + + def week_total_rule(model, w): + return ( + model.d_week_total[w] + model.a_min_week_total[w], + model.y_week_total[w], + model.d_week_total[w] + model.a_max_week_total[w], + ) + + # Constraint: Total weekly deliveries across all plants stay within tolerance band. + model.week_total_tolerance = pyo.Constraint(model.W, rule=week_total_rule) + + def month_total_rule(model, m): + return ( + model.d_month_total[m] + model.a_min_month_total[m], + model.y_month_total[m], + model.d_month_total[m] + model.a_max_month_total[m], + ) + + # Constraint: Total monthly deliveries across all plants stay within tolerance band. + model.month_total_tolerance = pyo.Constraint(model.M, rule=month_total_rule) + + model.a_min_V_day = pyo.Param(["N", "W", "ALL"], model.W, model.D, initialize=0.0, mutable=True) + model.a_max_V_day = pyo.Param(["N", "W", "ALL"], model.W, model.D, initialize=0.0, mutable=True) + model.a_min_V_week = pyo.Param(["N", "W", "ALL"], model.W, initialize=0.0, mutable=True) + model.a_max_V_week = pyo.Param(["N", "W", "ALL"], model.W, initialize=0.0, mutable=True) + model.a_min_V_month = pyo.Param(["N", "W", "ALL"], model.M, initialize=0.0, mutable=True) + model.a_max_V_month = pyo.Param(["N", "W", "ALL"], model.M, initialize=0.0, mutable=True) + + for w in model.W: + for d in model.D: + dN = model.dV_N[w, d] + dW = model.dV_W[w, d] + dA = dN + dW + model.a_min_V_day["N", w, d], model.a_max_V_day["N", w, d] = _tol_from_row( + ver_row_day_N, dN, step=step_size_tonnes + ) + model.a_min_V_day["W", w, d], model.a_max_V_day["W", w, d] = _tol_from_row( + ver_row_day_W, dW, step=step_size_tonnes + ) + model.a_min_V_day["ALL", w, d], model.a_max_V_day["ALL", w, d] = _tol_from_row( + ver_row_day_ALL, dA, step=step_size_tonnes + ) + + for w in model.W: + dN = sum(model.dV_N[w, dd] for dd in model.D) + dW = sum(model.dV_W[w, dd] for dd in model.D) + dA = dN + dW + a_min_n, a_max_n = _tol_from_row(ver_row_week_N, dN, step=step_size_tonnes) + a_min_w, a_max_w = _tol_from_row(ver_row_week_W, dW, step=step_size_tonnes) + a_min_a, a_max_a = _tol_from_row(ver_row_week_ALL, dA, step=step_size_tonnes) + model.a_min_V_week["N", w] = a_min_n - week_tolerance_debug + model.a_max_V_week["N", w] = a_max_n + week_tolerance_debug + model.a_min_V_week["W", w] = a_min_w - week_tolerance_debug + model.a_max_V_week["W", w] = a_max_w + week_tolerance_debug + model.a_min_V_week["ALL", w] = a_min_a - week_tolerance_debug + model.a_max_V_week["ALL", w] = a_max_a + week_tolerance_debug + + for m in model.M: + weeks_in_m = [w for w, mm in week_to_month.items() if mm == m and w in model.W] + dN = sum(model.dV_N[w, dd] for w in weeks_in_m for dd in model.D) + dW = sum(model.dV_W[w, dd] for w in weeks_in_m for dd in model.D) + dA = dN + dW + model.a_min_V_month["N", m], model.a_max_V_month["N", m] = _tol_from_row( + ver_row_month_N, dN, step=step_size_tonnes + ) + model.a_min_V_month["W", m], model.a_max_V_month["W", m] = _tol_from_row( + ver_row_month_W, dW, step=step_size_tonnes + ) + model.a_min_V_month["ALL", m], model.a_max_V_month["ALL", m] = _tol_from_row( + ver_row_month_ALL, dA, step=step_size_tonnes + ) + + def v_w_day_rule(m, w, d): + return ( + m.dV_W[w, d] + m.a_min_V_day["W", w, d], + sum(m.x[i, "V", w, d, s] for i in I_W for s in m.S), + m.dV_W[w, d] + m.a_max_V_day["W", w, d], + ) + + # Constraint: Veredlung day tolerance for Welzow coal. + model.v_w_day_tol = pyo.Constraint(model.W, model.D, rule=v_w_day_rule) + + def v_n_day_rule(m, w, d): + return ( + m.dV_N[w, d] + m.a_min_V_day["N", w, d], + sum(m.x[i, "V", w, d, s] for i in I_N for s in m.S), + m.dV_N[w, d] + m.a_max_V_day["N", w, d], + ) + + # Constraint: Veredlung day tolerance for Nochten coal. + model.v_n_day_tol = pyo.Constraint(model.W, model.D, rule=v_n_day_rule) + + def v_all_day_rule(m, w, d): + return ( + m.dV_W[w, d] + m.dV_N[w, d] + m.a_min_V_day["ALL", w, d], + sum(m.x[i, "V", w, d, s] for i in I_W + I_N for s in m.S), + m.dV_W[w, d] + m.dV_N[w, d] + m.a_max_V_day["ALL", w, d], + ) + + # Constraint: Veredlung day tolerance for combined coal types. + model.v_all_day_tol = pyo.Constraint(model.W, model.D, rule=v_all_day_rule) + + def v_w_week_rule(m, w): + return ( + sum(m.dV_W[w, d] for d in m.D) + m.a_min_V_week["W", w], + sum(m.x[i, "V", w, d, s] for i in I_W for d in m.D for s in m.S), + sum(m.dV_W[w, d] for d in m.D) + m.a_max_V_week["W", w], + ) + + # Constraint: Veredlung weekly tolerance for Welzow coal. + model.v_w_week_tol = pyo.Constraint(model.W, rule=v_w_week_rule) + + def v_n_week_rule(m, w): + return ( + sum(m.dV_N[w, d] for d in m.D) + m.a_min_V_week["N", w], + sum(m.x[i, "V", w, d, s] for i in I_N for d in m.D for s in m.S), + sum(m.dV_N[w, d] for d in m.D) + m.a_max_V_week["N", w], + ) + + # Constraint: Veredlung weekly tolerance for Nochten coal. + model.v_n_week_tol = pyo.Constraint(model.W, rule=v_n_week_rule) + + def v_all_week_rule(m, w): + return ( + sum(m.dV_W[w, d] + m.dV_N[w, d] for d in m.D) + m.a_min_V_week["ALL", w], + sum(m.x[i, "V", w, d, s] for i in I_W + I_N for d in m.D for s in m.S), + sum(m.dV_W[w, d] + m.dV_N[w, d] for d in m.D) + m.a_max_V_week["ALL", w], + ) + + # Constraint: Veredlung weekly tolerance for combined coal types. + model.v_all_week_tol = pyo.Constraint(model.W, rule=v_all_week_rule) + + def v_w_month_rule(m, mm): + weeks_in_m = [w for w, month in week_to_month.items() if month == mm and w in m.W] + return ( + sum(m.dV_W[w, d] for w in weeks_in_m for d in m.D) + m.a_min_V_month["W", mm], + sum(m.x[i, "V", w, d, s] for i in I_W for w in weeks_in_m for d in m.D for s in m.S), + sum(m.dV_W[w, d] for w in weeks_in_m for d in m.D) + m.a_max_V_month["W", mm], + ) + + # Constraint: Veredlung monthly tolerance for Welzow coal. + model.v_w_month_tol = pyo.Constraint(model.M, rule=v_w_month_rule) + + def v_n_month_rule(m, mm): + weeks_in_m = [w for w, month in week_to_month.items() if month == mm and w in m.W] + return ( + sum(m.dV_N[w, d] for w in weeks_in_m for d in m.D) + m.a_min_V_month["N", mm], + sum(m.x[i, "V", w, d, s] for i in I_N for w in weeks_in_m for d in m.D for s in m.S), + sum(m.dV_N[w, d] for w in weeks_in_m for d in m.D) + m.a_max_V_month["N", mm], + ) + + # Constraint: Veredlung monthly tolerance for Nochten coal. + model.v_n_month_tol = pyo.Constraint(model.M, rule=v_n_month_rule) + + def v_all_month_rule(m, mm): + weeks_in_m = [w for w, month in week_to_month.items() if month == mm and w in m.W] + return ( + sum(m.dV_W[w, d] + m.dV_N[w, d] for w in weeks_in_m for d in m.D) + m.a_min_V_month["ALL", mm], + sum(m.x[i, "V", w, d, s] for i in I_W + I_N for w in weeks_in_m for d in m.D for s in m.S), + sum(m.dV_W[w, d] + m.dV_N[w, d] for w in weeks_in_m for d in m.D) + m.a_max_V_month["ALL", mm], + ) + + # Constraint: Veredlung monthly tolerance for combined coal types. + model.v_all_month_tol = pyo.Constraint(model.M, rule=v_all_month_rule) + + j_map = { + "Jänschwalde": "J", + "Schwarze Pumpe": "SP", + "Boxberg Werk 3": "B3", + "Boxberg Werk 4": "B4", + } + + i_map = { + "Reichwalder-Kohle": "Reichwalde", + "Nochtener-Kohle": "Nochten", + "Welzower-Kohle": "Welzow", + } + + mix_df = tables["kohle_mix"] + + model.alpha_min = pyo.Param(model.I, model.J, initialize=0.0, mutable=True) + model.alpha_max = pyo.Param(model.I, model.J, initialize=1.0, mutable=True) + + for _, row in mix_df.iterrows(): + j = j_map.get(row["kraftwerk"]) + i = i_map.get(row["kohlesorte"]) + if (i in model.I) and (j in model.J): + model.alpha_min[i, j] = float(row["minimal"]) + model.alpha_max[i, j] = float(row["maximal"]) + + J_MIX = ["J", "SP", "B3", "B4", "V"] + + def x_day_rule(m, i, j, w, d): + return sum(m.x[i, j, w, d, s] for s in m.S) + + model.x_day = pyo.Expression(model.I, model.J, model.W, model.D, rule=x_day_rule) + + def y_day_rule(m, j, w, d): + return sum(m.x_day[i, j, w, d] for i in m.I) + + model.y_day = pyo.Expression(model.J, model.W, model.D, rule=y_day_rule) + + def mix_lower_rule(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + return pyo.Constraint.Skip + return m.alpha_min[i, j] * m.y_day[j, w, d] <= m.x_day[i, j, w, d] + + def mix_upper_rule(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + return pyo.Constraint.Skip + return m.x_day[i, j, w, d] <= m.alpha_max[i, j] * m.y_day[j, w, d] + + # Constraint: Mix lower/upper bound per source and plant (hard), per day. + model.mix_lower = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_lower_rule) + model.mix_upper = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_upper_rule) + + model.alpha_target_low = pyo.Param(model.I, model.J, initialize=0.0, mutable=True) + model.alpha_target_high = pyo.Param(model.I, model.J, initialize=1.0, mutable=True) + model.alpha_target_mid = pyo.Param(model.I, model.J, initialize=0.5, mutable=True) + + for _, row in mix_df.iterrows(): + j = j_map.get(row["kraftwerk"]) + i = i_map.get(row["kohlesorte"]) + if (i in model.I) and (j in model.J): + low = float(row["ziel_low"]) + high = float(row["ziel_high"]) + model.alpha_target_low[i, j] = low + model.alpha_target_high[i, j] = high + model.alpha_target_mid[i, j] = 0.5 * (low + high) + + model.mix_target_low_dev = pyo.Var(model.I, model.J, model.W, model.D, domain=pyo.NonNegativeReals) + model.mix_target_high_dev = pyo.Var(model.I, model.J, model.W, model.D, domain=pyo.NonNegativeReals) + model.mix_target_mid_dev = pyo.Var(model.I, model.J, model.W, model.D, domain=pyo.NonNegativeReals) + + def mix_target_low_rule(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + return pyo.Constraint.Skip + return m.x_day[i, j, w, d] + m.mix_target_low_dev[i, j, w, d] >= m.alpha_target_low[i, j] * m.y_day[j, w, d] + + def mix_target_high_rule(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + return pyo.Constraint.Skip + return m.x_day[i, j, w, d] <= m.alpha_target_high[i, j] * m.y_day[j, w, d] + m.mix_target_high_dev[i, j, w, d] + + # Constraint: Mix target lower/upper band (soft via deviation variable). + model.mix_target_low = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_target_low_rule) + model.mix_target_high = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_target_high_rule) + + def mix_target_mid_hi(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + return pyo.Constraint.Skip + return m.x_day[i, j, w, d] - m.alpha_target_mid[i, j] * m.y_day[j, w, d] <= m.mix_target_mid_dev[i, j, w, d] + + def mix_target_mid_lo(m, i, j, w, d): + if bunker_df is not None and (j not in J_MIX or j not in m.J_BUNKER): + return pyo.Constraint.Skip + return m.alpha_target_mid[i, j] * m.y_day[j, w, d] - m.x_day[i, j, w, d] <= m.mix_target_mid_dev[i, j, w, d] + + model.mix_target_mid_hi = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_target_mid_hi) + model.mix_target_mid_lo = pyo.Constraint(model.I, model.J, model.W, model.D, rule=mix_target_mid_lo) + + if bunker_df is not None: + vorfahr_df = tables.get("bunker_vorfahrfenster") + if vorfahr_df is not None and not vorfahr_df.empty: + vorfahrfenster_tage = int(vorfahr_df.iloc[0]["vorfahrfenster_tage"]) + else: + vorfahrfenster_tage = 0 + + date_to_wd = {pd.to_datetime(v).normalize(): k for k, v in wd_to_date.items()} + + model.bunker_target_short = pyo.Var(model.J_BUNKER, model.W, model.D, domain=pyo.NonNegativeReals) + + model.bunker_target_con = pyo.ConstraintList() + model.bunker_max_con = pyo.ConstraintList() + + for j in model.J_BUNKER: + for w in model.W: + for d in model.D: + model.bunker_max_con.add(model.bunker_total[j, w, d] <= model.bunker_max[j]) + + # B3 bunker mix targets (soft) to discourage Welzow and favor RW/NO in bunker. + model.b3_bunker_mix_low_dev = pyo.Var(model.I, model.W, model.D, domain=pyo.NonNegativeReals) + model.b3_bunker_mix_high_dev = pyo.Var(model.I, model.W, model.D, domain=pyo.NonNegativeReals) + + def b3_bunker_mix_low_rule(m, i, w, d): + if "B3" not in m.J_BUNKER: + return pyo.Constraint.Skip + return ( + m.bunker[i, "B3", w, d] + m.b3_bunker_mix_low_dev[i, w, d] + >= m.alpha_target_low[i, "B3"] * m.bunker_total["B3", w, d] + ) + + def b3_bunker_mix_high_rule(m, i, w, d): + if "B3" not in m.J_BUNKER: + return pyo.Constraint.Skip + return ( + m.bunker[i, "B3", w, d] + <= m.alpha_target_high[i, "B3"] * m.bunker_total["B3", w, d] + m.b3_bunker_mix_high_dev[i, w, d] + ) + + model.b3_bunker_mix_low = pyo.Constraint(model.I, model.W, model.D, rule=b3_bunker_mix_low_rule) + model.b3_bunker_mix_high = pyo.Constraint(model.I, model.W, model.D, rule=b3_bunker_mix_high_rule) + + # Zielbestand vor Montag (Sonntag) als soft constraint. + for date, (w, d) in date_to_wd.items(): + if date.strftime("%a") != "Mon": + continue + prev_date = date - pd.Timedelta(days=1) + if prev_date not in date_to_wd: + continue + w_prev, d_prev = date_to_wd[prev_date] + for j in model.J_BUNKER: + model.bunker_target_con.add( + model.bunker_target_short[j, w_prev, d_prev] + >= model.bunker_target[j] - model.bunker_total[j, w_prev, d_prev] + ) + + ih_map = {"welzow": ih_dates_welzow, "boxberg": ih_dates_boxberg} + + for src, ih_dates in ih_map.items(): + for ih_date in ih_dates: + target_date = ih_date - pd.Timedelta(days=vorfahrfenster_tage) + if target_date not in date_to_wd: + continue + w_prev, d_prev = date_to_wd[target_date] + for j in model.J_BUNKER: + if j not in J_MIX: + continue + model.bunker_target_con.add( + model.bunker_target_short[j, w_prev, d_prev] + >= model.bunker_target[j] - model.bunker_total[j, w_prev, d_prev] + ) + + model.lambda_J = pyo.Param(initialize=100_000, mutable=True, within=pyo.NonNegativeReals) + + def welzow_to_J(m): + return sum( + m.x["Welzow", "J", w, d, s] + for w in m.W + for d in m.D + for s in m.S + if ("Welzow", "J", w, d, s) in m.x + ) + + model.bonus_welzow_J = pyo.Expression(rule=welzow_to_J) + + def demand_rule(m, j, w, d): + return m.dV_N[w, d] + m.dV_W[w, d] if j == "V" else m.d[j, w, d] + + model.demand = pyo.Expression(model.J, model.W, model.D, rule=demand_rule) + model.dev = pyo.Expression( + model.J, model.W, model.D, rule=lambda m, j, w, d: m.y_delivery[j, w, d] - m.demand[j, w, d] + ) + + J_dev = list(J_power) + model.dev_pos = pyo.Var(J_dev, model.W, model.D, domain=pyo.NonNegativeReals) + model.dev_neg = pyo.Var(J_dev, model.W, model.D, domain=pyo.NonNegativeReals) + + def dev_balance_rule(model, j, w, d): + return model.dev_pos[j, w, d] - model.dev_neg[j, w, d] == model.dev[j, w, d] + + # Constraint: Split demand deviation into positive/negative components. + model.dev_balance = pyo.Constraint(J_dev, model.W, model.D, rule=dev_balance_rule) + + model.z_max = pyo.Var(J_power, model.W, model.D, domain=pyo.Binary) + + def max_reached_rule(m, j, w, d): + return m.dev_pos[j, w, d] <= m.a_max_day[j, w, d] * m.z_max[j, w, d] + + # Constraint: Flag when daily deviation hits the maximum tolerance. + model.max_reached = pyo.Constraint(J_power, model.W, model.D, rule=max_reached_rule) + + D_list = list(model.D.ordered_data()) + + def no_three_in_a_row_rule(m, j, w, idx): + d1, d2, d3 = D_list[idx : idx + 3] + return m.z_max[j, w, d1] + m.z_max[j, w, d2] + m.z_max[j, w, d3] <= 3 + + # Constraint: Prevent three consecutive days with max deviation. + model.no_three_in_a_row = pyo.Constraint(J_power, model.W, range(len(D_list) - 2), rule=no_three_in_a_row_rule) + + def y_V_nochten(m, w, d): + return sum(m.x["Nochten", "V", w, d, s] for s in m.S) + + def y_V_welzow(m, w, d): + return sum(m.x["Welzow", "V", w, d, s] for s in m.S) + + model.devV_N = pyo.Expression(model.W, model.D, rule=lambda m, w, d: y_V_nochten(m, w, d) - m.dV_N[w, d]) + model.devV_W = pyo.Expression(model.W, model.D, rule=lambda m, w, d: y_V_welzow(m, w, d) - m.dV_W[w, d]) + + model.devV_N_pos = pyo.Var(model.W, model.D, domain=pyo.NonNegativeReals) + model.devV_N_neg = pyo.Var(model.W, model.D, domain=pyo.NonNegativeReals) + model.devV_W_pos = pyo.Var(model.W, model.D, domain=pyo.NonNegativeReals) + model.devV_W_neg = pyo.Var(model.W, model.D, domain=pyo.NonNegativeReals) + + # Constraint: Split Nochten Veredlung deviation into positive/negative parts. + model.devV_N_balance = pyo.Constraint( + model.W, model.D, rule=lambda m, w, d: m.devV_N_pos[w, d] - m.devV_N_neg[w, d] == m.devV_N[w, d] + ) + # Constraint: Split Welzow Veredlung deviation into positive/negative parts. + model.devV_W_balance = pyo.Constraint( + model.W, model.D, rule=lambda m, w, d: m.devV_W_pos[w, d] - m.devV_W_neg[w, d] == m.devV_W[w, d] + ) + + model.lambda_glatt = pyo.Param(initialize=1e5, mutable=True, within=pyo.NonNegativeReals) + + model.m_ijwd = pyo.Expression( + model.I, + model.J, + model.W, + model.D, + rule=lambda m, i, j, w, d: (1 / len(m.S)) + * sum(m.x[i, j, w, d, s] for s in m.S if (i, j, w, d, s) in m.x), + ) + + model.t_glatt = pyo.Var(model.I, model.J, model.W, model.D, model.S, domain=pyo.NonNegativeReals) + + def glatt_hi(m, i, j, w, d, s): + return m.t_glatt[i, j, w, d, s] >= m.x[i, j, w, d, s] - m.m_ijwd[i, j, w, d] + + def glatt_lo(m, i, j, w, d, s): + return m.t_glatt[i, j, w, d, s] >= -(m.x[i, j, w, d, s] - m.m_ijwd[i, j, w, d]) + + # Constraint: Upper absolute deviation from average shift flow (smoothness). + model.glatt_hi = pyo.Constraint(model.I, model.J, model.W, model.D, model.S, rule=glatt_hi) + # Constraint: Lower absolute deviation from average shift flow (smoothness). + model.glatt_lo = pyo.Constraint(model.I, model.J, model.W, model.D, model.S, rule=glatt_lo) + + lambda_dev = 100_000 + lambda_v = 100_000 + lambda_shift = 100_000_000 + lambda_mix = 4_000_000_000 + lambda_mix_mid = 4_000_000 + lambda_b3_bunker_mix = 100_000_000 + lambda_shift_balance = 5_000_000 + lambda_bunker_target = 50_000_000 + SHIFT_TOL = step_size_tonnes + SHIFT_PAIRS = [("F", "S"), ("S", "N"), ("F", "N")] + model.SHIFT_PAIRS = pyo.Set(initialize=SHIFT_PAIRS, dimen=2) + + model.shift_dev_soft = pyo.Var(model.I, model.J, model.W, model.D, model.SHIFT_PAIRS, domain=pyo.NonNegativeReals) + model.shift_excess = pyo.Var(model.I, model.J, model.W, model.D, model.SHIFT_PAIRS, domain=pyo.NonNegativeReals) + + # ConstraintList: Enforce absolute shift differences and excess over tolerance. + model.shift_dev_abs = pyo.ConstraintList() + model.shift_excess_def = pyo.ConstraintList() + for i in model.I: + for j in model.J: + for w in model.W: + for d in model.D: + for s1, s2 in SHIFT_PAIRS: + if (i, j, w, d, s1) in model.x and (i, j, w, d, s2) in model.x: + model.shift_dev_abs.add( + model.shift_dev_soft[i, j, w, d, s1, s2] + >= model.x[i, j, w, d, s1] - model.x[i, j, w, d, s2] + ) + model.shift_dev_abs.add( + model.shift_dev_soft[i, j, w, d, s1, s2] + >= model.x[i, j, w, d, s2] - model.x[i, j, w, d, s1] + ) + model.shift_excess_def.add( + model.shift_excess[i, j, w, d, s1, s2] + >= model.shift_dev_soft[i, j, w, d, s1, s2] - SHIFT_TOL + ) + + def objective_rule(model): + deviation_penalty = ( + lambda_dev + * sum(model.dev_pos[j, w, d] + model.dev_neg[j, w, d] for j in J_dev for w in model.W for d in model.D) + + lambda_v + * sum( + model.devV_N_pos[w, d] + model.devV_N_neg[w, d] + model.devV_W_pos[w, d] + model.devV_W_neg[w, d] + for w in model.W + for d in model.D + ) + ) + + smoothness_penalty = lambda_shift * sum( + model.shift_excess[i, j, w, d, s1, s2] + for i in model.I + for j in model.J + for w in model.W + for d in model.D + for s1, s2 in model.SHIFT_PAIRS + if (i, j, w, d, s1) in model.x and (i, j, w, d, s2) in model.x + ) + + shift_balance_penalty = lambda_shift_balance * sum( + model.shift_balance_dev[j, w, d, s] for j in model.J for w in model.W for d in model.D for s in model.S + ) + + mix_target_penalty = lambda_mix * sum( + model.mix_target_low_dev[i, j, w, d] + model.mix_target_high_dev[i, j, w, d] + for i in model.I + for j in model.J + for w in model.W + for d in model.D + ) + + mix_target_mid_penalty = lambda_mix_mid * sum( + model.mix_target_mid_dev[i, j, w, d] for i in model.I for j in model.J for w in model.W for d in model.D + ) + + bunker_penalty = 0 + if bunker_df is not None: + bunker_penalty = lambda_bunker_target * sum( + model.bunker_target_short[j, w, d] for j in model.J_BUNKER for w in model.W for d in model.D + ) + b3_bunker_mix_penalty = lambda_b3_bunker_mix * sum( + model.b3_bunker_mix_low_dev[i, w, d] + model.b3_bunker_mix_high_dev[i, w, d] + for i in model.I + for w in model.W + for d in model.D + ) + else: + b3_bunker_mix_penalty = 0 + + return ( + deviation_penalty + + smoothness_penalty + + shift_balance_penalty + + mix_target_penalty + + mix_target_mid_penalty + + bunker_penalty + + b3_bunker_mix_penalty + ) + + # Objective: penalize deviations/smoothness, apply route bonuses and penalties. + model.obj = pyo.Objective(rule=objective_rule, sense=pyo.minimize) + + def forbid_reichwalde_V_rule(m, w, d, s): + return m.x["Reichwalde", "V", w, d, s] == 0 + + # Constraint: Disallow Reichwalde coal to Veredlung. + model.forbid_reichwalde_V = pyo.Constraint(model.W, model.D, model.S, rule=forbid_reichwalde_V_rule) + + def forbid_welzow_B4_rule(m, w, d, s): + return m.x["Welzow", "B4", w, d, s] == 0 + + # Constraint: Disallow Welzow coal to Boxberg Werk 4. + model.forbid_welzow_B4 = pyo.Constraint(model.W, model.D, model.S, rule=forbid_welzow_B4_rule) + + cap_df = tables["foerderkapaz"] + cap_month = cap_df[cap_df["zeitraum"] == "pro Monat"] + + i_map = { + "Reichwalde (RW)": "Reichwalde", + "Nochten (NO)": "Nochten", + "Welzow-Süd": "Welzow", + } + + model.F_max_month = pyo.Param(model.I, initialize=1e12, mutable=True) + for _, row in cap_month.iterrows(): + i = i_map.get(row["tagebau"]) + if i in model.I: + model.F_max_month[i] = float(row["maximal"]) + + def F_month_rule(m, i, mm): + weeks_in_m = [w for w, mth in week_to_month.items() if mth == mm and w in m.W] + return sum(m.x[i, j, w, d, s] for j in m.J for w in weeks_in_m for d in m.D for s in m.S) + + model.F_month = pyo.Expression(model.I, model.M, rule=F_month_rule) + + # Constraint: Monthly production cap per mine. + model.cap_month = pyo.Constraint(model.I, model.M, rule=lambda m, i, mm: m.F_month[i, mm] <= m.F_max_month[i]) + + # Constraint: Joint monthly cap for Reichwalde + Nochten. + model.cap_month_RWNO = pyo.Constraint( + model.M, rule=lambda m, mm: m.F_month["Reichwalde", mm] + m.F_month["Nochten", mm] <= 3_000_000.0 + ) + + verladung_cap = tables["verladungskap"] + + def cap_lookup(verladung_label: str, zeitraum_label: str) -> float: + row = verladung_cap[ + (verladung_cap["verladung"] == verladung_label) & (verladung_cap["zeitraum"] == zeitraum_label) + ] + if row.empty: + raise ValueError(f"Keine Verladungskap für {verladung_label} / {zeitraum_label}") + return float(row.iloc[0]["maximal"]) + + BOXBERG_SOURCES = ["Reichwalde", "Nochten"] + BOXBERG_DEST = ["J", "SP", "B3", "V"] + BOXBERG_SHIFT_CAP = cap_lookup("Boxberg (RW+NO)", "pro Schicht") + BOXBERG_DAY_CAP = cap_lookup("Boxberg (RW+NO)", "pro Tag") + + WELZOW_SOURCES = ["Welzow"] + WELZOW_DEST = ["J", "SP", "B3", "V"] + WELZOW_SHIFT_CAP = cap_lookup("Welzow-Süd", "pro Schicht") + WELZOW_DAY_CAP = cap_lookup("Welzow-Süd", "pro Tag") + + def boxberg_shift_cap_rule(m, w, d, s): + return sum(m.x[i, j, w, d, s] for i in BOXBERG_SOURCES if i in m.I for j in BOXBERG_DEST if j in m.J) <= BOXBERG_SHIFT_CAP + + # Constraint: Boxberg loading capacity per shift. + model.cap_boxberg_shift = pyo.Constraint(model.W, model.D, model.S, rule=boxberg_shift_cap_rule) + + def boxberg_day_cap_rule(m, w, d): + return sum( + m.x[i, j, w, d, s] + for i in BOXBERG_SOURCES + if i in m.I + for j in BOXBERG_DEST + if j in m.J + for s in m.S + ) <= BOXBERG_DAY_CAP + + # Constraint: Boxberg loading capacity per day. + model.cap_boxberg_day = pyo.Constraint(model.W, model.D, rule=boxberg_day_cap_rule) + + def welzow_shift_cap_rule(m, w, d, s): + return sum(m.x[i, j, w, d, s] for i in WELZOW_SOURCES if i in m.I for j in WELZOW_DEST if j in m.J) <= WELZOW_SHIFT_CAP + + # Constraint: Welzow loading capacity per shift. + model.cap_welzow_shift = pyo.Constraint(model.W, model.D, model.S, rule=welzow_shift_cap_rule) + + def welzow_day_cap_rule(m, w, d): + return sum( + m.x[i, j, w, d, s] + for i in WELZOW_SOURCES + if i in m.I + for j in WELZOW_DEST + if j in m.J + for s in m.S + ) <= WELZOW_DAY_CAP + + # Constraint: Welzow loading capacity per day. + model.cap_welzow_day = pyo.Constraint(model.W, model.D, rule=welzow_day_cap_rule) + + avail = tables.get("Verfuegbarkeiten") + if avail is not None: + day_map = {"Sat": "Sa", "Sun": "So", "Mon": "Mo", "Tue": "Di", "Wed": "Mi", "Thu": "Do", "Fri": "Fr"} + schicht_order = ["F", "S", "N"] + + model.cap_Welzow = pyo.Param( + model.W, model.D, model.S, initialize=math.nan, mutable=True, within=pyo.Any + ) + model.cap_RW_N = pyo.Param( + model.W, model.D, model.S, initialize=math.nan, mutable=True, within=pyo.Any + ) + + for _, row in avail.iterrows(): + datum = row["datum"] + w = _mining_week_from_date(datum) + d = day_map[datum.strftime("%a")] + + caps_w = [row["Welzow_Sued_S1_t"], row["Welzow_Sued_S2_t"], row["Welzow_Sued_S3_t"]] + caps_rn = [row["Boxberg_NO_RW_S1_t"], row["Boxberg_NO_RW_S2_t"], row["Boxberg_NO_RW_S3_t"]] + + if all((not math.isnan(v)) and v == 0 for v in caps_w): + ih_dates_welzow.add(pd.to_datetime(datum).normalize()) + if all((not math.isnan(v)) and v == 0 for v in caps_rn): + ih_dates_boxberg.add(pd.to_datetime(datum).normalize()) + + for idx, s in enumerate(schicht_order): + val_w = caps_w[idx] + if not math.isnan(val_w): + model.cap_Welzow[w, d, s] = val_w + val_rn = caps_rn[idx] + if not math.isnan(val_rn): + model.cap_RW_N[w, d, s] = val_rn + + def cap_welzow_rule(m, w, d, s): + v = value(m.cap_Welzow[w, d, s]) + if math.isnan(v): + return pyo.Constraint.Skip + return ( + sum( + m.x[i, j, w, d, s] + for i in ["Welzow"] + for j in m.J + if (i, j, w, d, s) in m.x.index_set() + ) + <= v + ) + + # Constraint: Dynamic availability cap for Welzow (per shift). + model.cap_welzow_con = pyo.Constraint(model.W, model.D, model.S, rule=cap_welzow_rule) + + def cap_rw_n_rule(m, w, d, s): + v = value(m.cap_RW_N[w, d, s]) + if math.isnan(v): + return pyo.Constraint.Skip + return ( + sum( + m.x[i, j, w, d, s] + for i in ["Reichwalde", "Nochten"] + for j in m.J + if (i, j, w, d, s) in m.x.index_set() + ) + <= v + ) + + # Constraint: Dynamic availability cap for Reichwalde + Nochten (per shift). + model.cap_rw_n_con = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_rule) + + zug = tables.get("zugdurchlass") + if zug is not None: + df = zug.copy() + + def cap_from_table(start_exact: str, ziel_exact: str) -> float: + s = df["start"].astype(str).str.strip() + z = df["ziel"].astype(str).str.strip() + mask = (s == start_exact) & (z == ziel_exact) + row = df[mask] + if row.empty: + raise ValueError(f"cap not found for start={start_exact!r}, ziel={ziel_exact!r}") + return float(row.iloc[0]["maximal"]) + + CAP_RW_N_TO_J = cap_from_table("Verladung Boxberg (KLP)", "KW Jänschwalde") + CAP_RW_N_TO_SP = cap_from_table("Verladung Boxberg (KLP)", "KW Schwarze Pumpe") + CAP_RW_N_TO_V = cap_from_table("Verladung Boxberg (KLP)", "Veredlung ISP") + CAP_RW_N_TO_B3 = cap_from_table("Verladung Boxberg (KLP)", "KW Boxberg Werk 3") + + CAP_W_TO_J = cap_from_table("Verladung Welzow-Süd (KUP)", "KW Jänschwalde") + CAP_W_TO_SP = cap_from_table("Verladung Welzow-Süd (KUP)", "KW Schwarze Pumpe") + CAP_W_TO_V = cap_from_table("Verladung Welzow-Süd (KUP)", "Veredlung ISP") + CAP_W_TO_B3 = cap_from_table("Verladung Welzow-Süd (KUP)", "KW Boxberg Werk 3") + + CAP_RW_N_TO_SP_V = cap_from_table("Verladung Boxberg (KLP)", "KW SP + V ISP") + CAP_RW_N_TO_ALL = cap_from_table("Verladung Boxberg (KLP)", "KW JW + KW SP + V ISP + KW B3") + + CAP_W_TO_SP_V = cap_from_table("Verladung Welzow-Süd (KUP)", "KW SP + V ISP") + CAP_W_TO_SP_V_B3 = cap_from_table("Verladung Welzow-Süd (KUP)", "KW SP + V ISP + KW B3") + + CAP_RW_N_W_TO_J = cap_from_table("KUP + KLP", "KW Jänschwalde") + CAP_RW_N_W_TO_B3 = cap_from_table("KUP + KLP", "KW Boxberg Werk 3") + + CAP_RW_N_J_AND_W_B3 = cap_from_table("KLP zum KW JW", "KUP zum KW B3") + CAP_RW_N_JSPV_AND_W_B3 = cap_from_table("KLP zum KW JW + KW SP + V ISP", "KUP zum KW B3") + + def cap_rw_n_to_j_rule(m, w, d, s): + return m.x["Reichwalde", "J", w, d, s] + m.x["Nochten", "J", w, d, s] <= CAP_RW_N_TO_J + + # Constraint: KLP capacity from RW+NO to J. + model.cap_rw_n_to_j = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_j_rule) + + def cap_rw_n_to_sp_rule(m, w, d, s): + return m.x["Reichwalde", "SP", w, d, s] + m.x["Nochten", "SP", w, d, s] <= CAP_RW_N_TO_SP + + # Constraint: KLP capacity from RW+NO to SP. + model.cap_rw_n_to_sp = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_sp_rule) + + def cap_rw_n_to_v_rule(m, w, d, s): + return m.x["Reichwalde", "V", w, d, s] + m.x["Nochten", "V", w, d, s] <= CAP_RW_N_TO_V + + # Constraint: KLP capacity from RW+NO to Veredlung. + model.cap_rw_n_to_v = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_v_rule) + + def cap_rw_n_to_bw3_rule(m, w, d, s): + return m.x["Reichwalde", "B3", w, d, s] + m.x["Nochten", "B3", w, d, s] <= CAP_RW_N_TO_B3 + + # Constraint: KLP capacity from RW+NO to B3. + model.cap_rw_n_to_bw3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_bw3_rule) + + def cap_w_to_j_rule(m, w, d, s): + return m.x["Welzow", "J", w, d, s] <= CAP_W_TO_J + + # Constraint: KUP capacity from Welzow to J. + model.cap_w_to_j = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_j_rule) + + model.k_w_to_j_steps = pyo.Var(model.W, model.D, model.S, domain=pyo.NonNegativeIntegers) + + def welzow_j_multiple_rule(m, w, d, s): + return m.x["Welzow", "J", w, d, s] == (2 * m.step_size_tonnes) * m.k_w_to_j_steps[w, d, s] + + # Constraint: Enforce 2kt multiples for Welzow->J flows. + model.welzow_j_multiple_2kt = pyo.Constraint(model.W, model.D, model.S, rule=welzow_j_multiple_rule) + + def cap_w_to_sp_rule(m, w, d, s): + return m.x["Welzow", "SP", w, d, s] <= CAP_W_TO_SP + + # Constraint: KUP capacity from Welzow to SP. + model.cap_w_to_sp = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_sp_rule) + + def cap_w_to_v_rule(m, w, d, s): + return m.x["Welzow", "V", w, d, s] <= CAP_W_TO_V + + # Constraint: KUP capacity from Welzow to Veredlung. + model.cap_w_to_v = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_v_rule) + + def cap_w_to_bw3_rule(m, w, d, s): + return m.x["Welzow", "B3", w, d, s] <= CAP_W_TO_B3 + + # Constraint: KUP capacity from Welzow to B3. + model.cap_w_to_bw3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_bw3_rule) + + def cap_rw_n_to_sp_v_rule(m, w, d, s): + return ( + m.x["Reichwalde", "SP", w, d, s] + + m.x["Nochten", "SP", w, d, s] + + m.x["Reichwalde", "V", w, d, s] + + m.x["Nochten", "V", w, d, s] + <= CAP_RW_N_TO_SP_V + ) + + # Constraint: KLP combined capacity to SP + Veredlung. + model.cap_rw_n_to_sp_v = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_sp_v_rule) + + def cap_rw_n_to_j_sp_v_b3_rule(m, w, d, s): + return ( + m.x["Reichwalde", "J", w, d, s] + + m.x["Reichwalde", "SP", w, d, s] + + m.x["Reichwalde", "V", w, d, s] + + m.x["Reichwalde", "B3", w, d, s] + + m.x["Nochten", "J", w, d, s] + + m.x["Nochten", "SP", w, d, s] + + m.x["Nochten", "V", w, d, s] + + m.x["Nochten", "B3", w, d, s] + <= CAP_RW_N_TO_ALL + ) + + # Constraint: KLP combined capacity to J + SP + V + B3. + model.cap_rw_n_to_j_sp_v_b3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_to_j_sp_v_b3_rule) + + def cap_w_to_sp_v_rule(m, w, d, s): + return m.x["Welzow", "SP", w, d, s] + m.x["Welzow", "V", w, d, s] <= CAP_W_TO_SP_V + + # Constraint: KUP combined capacity to SP + V. + model.cap_w_to_sp_v = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_sp_v_rule) + + def cap_w_to_sp_v_b3_rule(m, w, d, s): + return ( + m.x["Welzow", "SP", w, d, s] + + m.x["Welzow", "V", w, d, s] + + m.x["Welzow", "B3", w, d, s] + <= CAP_W_TO_SP_V_B3 + ) + + # Constraint: KUP combined capacity to SP + V + B3. + model.cap_w_to_sp_v_b3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_w_to_sp_v_b3_rule) + + def cap_rw_n_w_to_j_rule(m, w, d, s): + return ( + m.x["Reichwalde", "J", w, d, s] + + m.x["Nochten", "J", w, d, s] + + m.x["Welzow", "J", w, d, s] + <= CAP_RW_N_W_TO_J + ) + + # Constraint: Combined KLP/KUP capacity to J. + model.cap_rw_n_w_to_j = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_w_to_j_rule) + + kvb_nord = tables.get("zugdurchlass_kvb_nord") + print(kvb_nord) + if kvb_nord is not None: + day_map = {"Sat": "Sa", "Sun": "So", "Mon": "Mo", "Tue": "Di", "Wed": "Mi", "Thu": "Do", "Fri": "Fr"} + schicht_order = ["F", "S", "N"] + + model.cap_KVB_Nord = pyo.Param( + model.W, + model.D, + model.S, + initialize=math.nan, + mutable=True, + within=pyo.Any, + ) + + for _, row in kvb_nord.iterrows(): + datum = row["datum"] + w = _mining_week_from_date(datum) + d = day_map[datum.strftime("%a")] + + caps = [row["KVB_Nord_S1_t"], row["KVB_Nord_S2_t"], row["KVB_Nord_S3_t"]] + for idx, s in enumerate(schicht_order): + val = caps[idx] + if not math.isnan(val): + model.cap_KVB_Nord[w, d, s] = val + + def cap_kvb_nord_rule(m, w, d, s): + v = value(m.cap_KVB_Nord[w, d, s]) + if math.isnan(v): + return pyo.Constraint.Skip + return ( + m.x["Reichwalde", "J", w, d, s] + + m.x["Nochten", "J", w, d, s] + + m.x["Welzow", "J", w, d, s] + <= v + ) + + # Constraint: KVB Nord short-term capacity to J (per day/shift if provided). + model.cap_kvb_nord = pyo.Constraint(model.W, model.D, model.S, rule=cap_kvb_nord_rule) + + def cap_rw_n_w_to_bw3_rule(m, w, d, s): + return ( + m.x["Reichwalde", "B3", w, d, s] + + m.x["Nochten", "B3", w, d, s] + + m.x["Welzow", "B3", w, d, s] + <= CAP_RW_N_W_TO_B3 + ) + + # Constraint: Combined KLP/KUP capacity to B3. + model.cap_rw_n_w_to_bw3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_w_to_bw3_rule) + + def cap_rw_n_j_and_w_b3_rule(m, w, d, s): + return ( + m.x["Reichwalde", "J", w, d, s] + + m.x["Nochten", "J", w, d, s] + + m.x["Welzow", "B3", w, d, s] + <= CAP_RW_N_J_AND_W_B3 + ) + + # Constraint: Combined capacity (KLP->J, KUP->B3). + model.cap_rw_n_j_and_w_b3 = pyo.Constraint(model.W, model.D, model.S, rule=cap_rw_n_j_and_w_b3_rule) + + def cap_rw_n_to_j_sp_v_and_w_to_b3_rule(m, w, d, s): + return ( + m.x["Reichwalde", "J", w, d, s] + + m.x["Reichwalde", "SP", w, d, s] + + m.x["Reichwalde", "V", w, d, s] + + m.x["Nochten", "J", w, d, s] + + m.x["Nochten", "SP", w, d, s] + + m.x["Nochten", "V", w, d, s] + + m.x["Welzow", "B3", w, d, s] + <= CAP_RW_N_JSPV_AND_W_B3 + ) + + # Constraint: Combined capacity (KLP->J/SP/V, KUP->B3). + model.cap_rw_n_to_j_sp_v_and_w_to_b3 = pyo.Constraint( + model.W, model.D, model.S, rule=cap_rw_n_to_j_sp_v_and_w_to_b3_rule + ) + + model._wd_to_date = wd_to_date + model._week_to_month = week_to_month + + return model + + +def solve_model( + model: pyo.ConcreteModel, + solver_name: str = "gurobi", + time_limit: int = 600, + mip_gap: float = 0.05, + iis_path: Path = Path("results/iis.ilp"), +) -> pyo.SolverResults: + if solver_name == "highs": + solver = pyo.SolverFactory("highs") + solver.options.update( + { + "time_limit": time_limit, + "mip_rel_gap": mip_gap, + "threads": 8, + "output_flag": True, + } + ) + return solver.solve(model, tee=True) + + opt = pyo.SolverFactory(solver_name) + if solver_name == "gurobi": + opt.options["TimeLimit"] = time_limit + opt.options["MIPGap"] = mip_gap + opt.options["LogFile"] = "gurobi_first_opt_model.log" + elif solver_name == "scip": + # SCIP uses hierarchical parameter names. + opt.options["limits/time"] = time_limit + opt.options["limits/gap"] = mip_gap + results = opt.solve(model, tee=True, symbolic_solver_labels=True) + + if solver_name == "gurobi" and results.solver.termination_condition in { + TerminationCondition.infeasible, + TerminationCondition.infeasibleOrUnbounded, + }: + try: + iis_path.parent.mkdir(parents=True, exist_ok=True) + persistent = pyo.SolverFactory("gurobi_persistent") + persistent.set_instance(model, symbolic_solver_labels=True) + persistent.options["TimeLimit"] = time_limit + persistent.options["MIPGap"] = mip_gap + persistent.solve(tee=False) + persistent._solver_model.computeIIS() + persistent._solver_model.write(str(iis_path)) + except Exception as exc: + print(f"WARNING: Failed to write IIS file: {exc}") + + return results diff --git a/src/optimization/run_optimization.py b/src/optimization/run_optimization.py new file mode 100644 index 0000000..e15588d --- /dev/null +++ b/src/optimization/run_optimization.py @@ -0,0 +1,731 @@ +from __future__ import annotations + +import argparse +import sys +from pathlib import Path + +import pandas as pd +import pyomo.environ as pyo +from pyomo.environ import value + +SRC_ROOT = Path(__file__).resolve().parents[1] +if str(SRC_ROOT) not in sys.path: + sys.path.insert(0, str(SRC_ROOT)) + +from optimization.model_builder import build_model, load_tables, solve_model + + + + +def report_results(model: pyo.ConcreteModel, max_rows: int) -> None: + print("Objective value:", value(model.obj)) + print("Non-zero production decisions (k):") + printed = 0 + for i in model.I: + for j in model.J: + for w in model.W: + for d in model.D: + for s in model.S: + qty = value(model.k[i, j, w, d, s]) + if qty > 1e-6: + print(f" {i} -> {j} (W{w} D{d} S{s}): {qty:.0f}") + printed += 1 + if printed >= max_rows: + print(" ... output truncated ...") + return + + +def export_results(model: pyo.ConcreteModel, output_path: Path) -> None: + output_path.parent.mkdir(parents=True, exist_ok=True) + wd_to_date = getattr(model, "_wd_to_date", {}) + + def safe_value(var) -> float: + val = pyo.value(var, exception=False) + return float(val) if val is not None else 0.0 + + def autosize_worksheet(ws, df, index_cols=None, max_width=25): + if index_cols is None: + index_cols = list(df.index.names) + idx_names = list(index_cols) + col_widths = [max(10, max(len(str(n)) for n in idx_names if n is not None) + 2) if idx_names else 10] + for col_idx, col in enumerate(df.columns): + header = " / ".join([str(c) for c in col]) if isinstance(col, tuple) else str(col) + max_len = max(len(header), 8) + sample = df.iloc[:200, col_idx] + max_len = max(max_len, sample.astype(str).str.len().max()) + col_widths.append(min(max_width, int(max_len) + 2)) + return col_widths + + def adjust_widths_for_labels(df, widths, label_scale, index_scale=None): + adjusted = widths[:] + if index_scale is not None and adjusted: + adjusted[0] = max(10, int(adjusted[0] * index_scale)) + if hasattr(df.columns, "get_level_values"): + top_level = df.columns.get_level_values(0) + for idx, label in enumerate(top_level, start=1): + if label in label_scale: + adjusted[idx] = max(6, int(adjusted[idx] * label_scale[label])) + return adjusted + lieferungen_schicht = [] + for j in model.J: + for w in model.W: + for d in model.D: + for s in model.S: + if j == "V": + nachfrage = pyo.value(model.dV_N[w, d] + model.dV_W[w, d]) + else: + nachfrage = pyo.value(model.d[j, w, d]) + use_bunker_out = hasattr(model, "bunker_out") and j in getattr(model, "J_BUNKER", []) + delivery_sum = sum(safe_value(model.x[i, j, w, d, s]) for i in model.I) + out_sum = ( + sum(safe_value(model.bunker_out[i, j, w, d, s]) for i in model.I) + if use_bunker_out + else delivery_sum + ) + bunker_inflow = 0.0 + if use_bunker_out: + x_sum = sum(safe_value(model.x[i, j, w, d, s]) for i in model.I) + out_sum = sum(safe_value(model.bunker_out[i, j, w, d, s]) for i in model.I) + bunker_inflow = round(x_sum - out_sum, 2) + def flow_val(i_name: str) -> float: + return safe_value(model.x[i_name, j, w, d, s]) + lieferungen_schicht.append( + { + "kraftwerk": j, + "woche": w, + "tag": d, + "datum": wd_to_date.get((w, d)), + "schicht": s, + "nachfrage_tonnen": nachfrage, + "lieferung_tonnen": delivery_sum, + "lieferungsabweichung_tonnen": round(delivery_sum - nachfrage, 2), + "bunkerzufluss_tonnen": bunker_inflow, + "Nochten": flow_val("Nochten"), + "Reichwalde": flow_val("Reichwalde"), + "Welzow": flow_val("Welzow"), + } + ) + + order_k_pw = ["J", "SP", "B3", "B4"] + order_k_v = ["V"] + order_sources = ["Reichwalde", "Nochten", "Welzow"] + order_s = ["F", "S", "N"] + + df_raw = pd.DataFrame(lieferungen_schicht).copy() + df_raw["datum"] = pd.to_datetime(df_raw["datum"]) + + v_demand_map = { + (int(w), d): { + "welzow": float(pyo.value(model.dV_W[w, d])), + "nochten": float(pyo.value(model.dV_N[w, d])), + } + for w in model.W + for d in model.D + } + + df_raw["nachfrage_welzow"] = df_raw.get("nachfrage_welzow", pd.Series(index=df_raw.index)) + df_raw["nachfrage_nochten"] = df_raw.get("nachfrage_nochten", pd.Series(index=df_raw.index)) + df_raw["nachfrage_welzow"] = df_raw.apply( + lambda r: ( + r["nachfrage_welzow"] + if pd.notna(r["nachfrage_welzow"]) + else v_demand_map.get((int(r["woche"]), r["tag"]), {}).get("welzow", 0) + ), + axis=1, + ) + df_raw["nachfrage_nochten"] = df_raw.apply( + lambda r: ( + r["nachfrage_nochten"] + if pd.notna(r["nachfrage_nochten"]) + else v_demand_map.get((int(r["woche"]), r["tag"]), {}).get("nochten", 0) + ), + axis=1, + ) + + df = df_raw.rename(columns={"lieferungen_tonnen": "lieferung_tonnen"}).copy() + if "lieferung_tonnen" not in df.columns: + df["lieferung_tonnen"] = df[order_sources].sum(axis=1) + + present_sources = [c for c in order_sources if c in df.columns] + + df_src = ( + df.pivot_table( + index=["datum", "woche", "tag"], + columns=["kraftwerk", "schicht"], + values=present_sources, + aggfunc="sum", + ) + .fillna(0) + ) + df_src.columns = df_src.columns.reorder_levels([1, 0, 2]) + df_src = df_src.reindex( + columns=pd.MultiIndex.from_product([order_k_pw + order_k_v, present_sources, order_s]), + fill_value=0, + ) + + df_demand = ( + df.groupby(["datum", "woche", "tag", "kraftwerk"])["nachfrage_tonnen"] + .first() + .unstack("kraftwerk") + .reindex(columns=order_k_pw + order_k_v, fill_value=0) + .reindex(df_src.index, fill_value=0) + ) + + df_v_demand_split = ( + df[df["kraftwerk"] == "V"] + .groupby(["datum", "woche", "tag"])[["nachfrage_welzow", "nachfrage_nochten"]] + .first() + .reindex(df_src.index, fill_value=0) + ) + + totals_plain = df_src.T.groupby(level=0).sum().T + totals_plain = totals_plain.reindex(columns=order_k_pw + order_k_v, fill_value=0) + + day_diff_rows = [] + for w in model.W: + for d in model.D: + date = wd_to_date.get((w, d)) + if date is None: + continue + row = {"datum": date, "woche": w, "tag": d} + for j in order_k_pw + order_k_v: + if j == "V": + demand = pyo.value(model.dV_N[w, d] + model.dV_W[w, d]) + else: + demand = pyo.value(model.d[j, w, d]) + delivered = pyo.value(model.y_delivery[j, w, d]) + row[j] = delivered - demand + day_diff_rows.append(row) + + day_diff_plain = ( + pd.DataFrame(day_diff_rows) + .set_index(["datum", "woche", "tag"]) + .reindex(totals_plain.index, fill_value=0) + .reindex(columns=order_k_pw + order_k_v, fill_value=0) + ) + + totals = totals_plain.copy() + totals.columns = pd.MultiIndex.from_tuples([(k, "Gesamt", "") for k in totals.columns]) + + demand_cols = df_demand.copy() + demand_cols.columns = pd.MultiIndex.from_tuples([(k, "Nachfrage", "") for k in demand_cols.columns]) + + day_diff_cols = day_diff_plain.copy() + day_diff_cols.columns = pd.MultiIndex.from_tuples( + [(k, "Lieferungstagesabweichung", "") for k in day_diff_cols.columns] + ) + + v_demand_cols = df_v_demand_split.copy() + v_demand_cols.columns = pd.MultiIndex.from_tuples( + [ + ("V", "Nachfrage_Welzow", ""), + ("V", "Nachfrage_Nochtener", ""), + ] + ) + + has_bunker = hasattr(model, "bunker") + col_order = [] + for k in order_k_pw: + for src in present_sources: + for sch in order_s: + col_order.append((k, src, sch)) + col_order.append((k, "Gesamt", "")) + col_order.append((k, "Nachfrage", "")) + col_order.append((k, "Lieferungstagesabweichung", "")) + + col_order += [( "V", src, sch) for src in present_sources for sch in order_s] + col_order += [ + ("V", "Nachfrage_Welzow", ""), + ("V", "Nachfrage_Nochtener", ""), + ("V", "Gesamt", ""), + ("V", "Nachfrage", ""), + ("V", "Lieferungstagesabweichung", ""), + ] + + df_out = pd.concat([df_src, v_demand_cols, totals, demand_cols, day_diff_cols], axis=1) + + df_out = df_out.reindex( + columns=col_order, + fill_value=0, + ) + + weekday_order = ["Mo", "Di", "Mi", "Do", "Fr", "Sa", "So"] + idx = df_out.index + df_out = df_out.copy() + df_out.index = pd.MultiIndex.from_arrays( + [ + pd.to_datetime(idx.get_level_values("datum")), + idx.get_level_values("woche"), + pd.Categorical(idx.get_level_values("tag"), categories=weekday_order, ordered=True), + ], + names=["datum", "woche", "tag"], + ) + df_out = df_out.sort_index(level=["datum", "woche", "tag"]) + + df = df_out.copy() + df_out = df_out / 1000 + + bunker_sheet = None + if has_bunker: + bunker_rows = [] + for j in getattr(model, "J_BUNKER", []): + for w in model.W: + for d in model.D: + date = wd_to_date.get((w, d)) + if date is None: + continue + bunker_total = sum(safe_value(model.bunker[i, j, w, d]) for i in model.I) + bunker_rows.append( + { + "kraftwerk": j, + "woche": w, + "tag": d, + "datum": date, + "bunkerbestand_tonnen": bunker_total, + } + ) + + if bunker_rows: + bunker_df = pd.DataFrame(bunker_rows) + bunker_df["datum"] = pd.to_datetime(bunker_df["datum"]) + bunker_df = bunker_df.sort_values(["kraftwerk", "datum", "woche", "tag"]) + bunker_df["vortags_bunkerbestand_tonnen"] = bunker_df.groupby("kraftwerk")[ + "bunkerbestand_tonnen" + ].shift(1) + bunker_df["vortags_bunkerbestand_tonnen"] = bunker_df["vortags_bunkerbestand_tonnen"].fillna(0) + bunker_pivot = ( + bunker_df.pivot_table( + index=["datum", "woche", "tag"], + columns=["kraftwerk"], + values=["bunkerbestand_tonnen"], + aggfunc="first", + ) + .fillna(0) + ) + bunker_pivot = bunker_pivot.reindex(columns=order_k_pw + order_k_v, level=1, fill_value=0) + bunker_pivot.columns = pd.MultiIndex.from_tuples( + [(k, "Bunkerbestand", "") for k in bunker_pivot.columns.get_level_values(1)] + ) + bunker_prev_pivot = ( + bunker_df.pivot_table( + index=["datum", "woche", "tag"], + columns=["kraftwerk"], + values=["vortags_bunkerbestand_tonnen"], + aggfunc="first", + ) + .fillna(0) + ) + bunker_prev_pivot = bunker_prev_pivot.reindex(columns=order_k_pw + order_k_v, level=1, fill_value=0) + bunker_prev_pivot.columns = pd.MultiIndex.from_tuples( + [(k, "Vortags_Bunkerbestand", "") for k in bunker_prev_pivot.columns.get_level_values(1)] + ) + inflow_pivot = None + if "bunkerzufluss_tonnen" in df_raw.columns: + inflow_pivot = ( + df_raw.pivot_table( + index=["datum", "woche", "tag"], + columns=["kraftwerk"], + values=["bunkerzufluss_tonnen"], + aggfunc="sum", + ) + .fillna(0) + ) + inflow_pivot.columns = pd.MultiIndex.from_tuples( + [(k, "Bunkerzufluss", "") for k in inflow_pivot.columns.get_level_values(1)] + ) + + frames = [df] + frames.append(bunker_prev_pivot.reindex(df.index, fill_value=0)) + if inflow_pivot is not None: + frames.append(inflow_pivot.reindex(df.index, fill_value=0)) + frames.append(bunker_pivot.reindex(df.index, fill_value=0)) + bunker_sheet = pd.concat(frames, axis=1) + + col_order_bunker = [] + for k in order_k_pw: + for src in present_sources: + for sch in order_s: + col_order_bunker.append((k, src, sch)) + col_order_bunker.append((k, "Gesamt", "")) + col_order_bunker.append((k, "Nachfrage", "")) + col_order_bunker.append((k, "Lieferungstagesabweichung", "")) + if (k, "Vortags_Bunkerbestand", "") in bunker_sheet.columns: + col_order_bunker.append((k, "Vortags_Bunkerbestand", "")) + if (k, "Bunkerzufluss", "") in bunker_sheet.columns: + col_order_bunker.append((k, "Bunkerzufluss", "")) + if (k, "Bunkerbestand", "") in bunker_sheet.columns: + col_order_bunker.append((k, "Bunkerbestand", "")) + + col_order_bunker += [("V", src, sch) for src in present_sources for sch in order_s] + col_order_bunker += [ + ("V", "Nachfrage_Welzow", ""), + ("V", "Nachfrage_Nochtener", ""), + ("V", "Gesamt", ""), + ("V", "Nachfrage", ""), + ("V", "Lieferungstagesabweichung", ""), + ] + if ("V", "Vortags_Bunkerbestand", "") in bunker_sheet.columns: + col_order_bunker.append(("V", "Vortags_Bunkerbestand", "")) + if ("V", "Bunkerzufluss", "") in bunker_sheet.columns: + col_order_bunker.append(("V", "Bunkerzufluss", "")) + if ("V", "Bunkerbestand", "") in bunker_sheet.columns: + col_order_bunker.append(("V", "Bunkerbestand", "")) + + bunker_sheet = bunker_sheet.reindex(columns=col_order_bunker, fill_value=0) + + try: + import xlsxwriter # type: ignore + + excel_engine = "xlsxwriter" + except Exception: + excel_engine = "openpyxl" + + with pd.ExcelWriter(output_path, engine=excel_engine) as writer: + df.to_excel(writer, sheet_name="Sheet1") + if excel_engine == "xlsxwriter": + ws1 = writer.sheets["Sheet1"] + widths = autosize_worksheet(ws1, df) + widths = adjust_widths_for_labels( + df, + widths, + { + "Reichwalde": 0.5, + "Nochten": 0.5, + "Welzow": 0.5, + "Gesamt": 0.5, + "Nachfrage": 0.5, + "Bunkerbestand": 0.5, + "Bunkerzufluss": 0.5, + "Vortags_Bunkerbestand": 0.5, + }, + index_scale=1.2, + ) + for i, w in enumerate(widths): + ws1.set_column(i, i, w) + else: + ws1 = writer.sheets["Sheet1"] + widths = autosize_worksheet(ws1, df) + widths = adjust_widths_for_labels( + df, + widths, + { + "Reichwalde": 0.5, + "Nochten": 0.5, + "Welzow": 0.5, + "Gesamt": 0.5, + "Nachfrage": 0.5, + "Bunkerbestand": 0.5, + "Bunkerzufluss": 0.5, + "Vortags_Bunkerbestand": 0.5, + }, + index_scale=1.2, + ) + for i, w in enumerate(widths): + ws1.column_dimensions[chr(65 + i)].width = w + + if bunker_sheet is not None: + bunker_sheet.to_excel(writer, sheet_name="mit_Bunkerbestand") + if excel_engine == "xlsxwriter": + workbook = writer.book + worksheet = writer.sheets["mit_Bunkerbestand"] + widths = autosize_worksheet(worksheet, bunker_sheet) + widths = adjust_widths_for_labels( + bunker_sheet, + widths, + { + "Reichwalde": 0.5, + "Nochten": 0.5, + "Welzow": 0.5, + "Gesamt": 0.5, + "Nachfrage": 0.5, + "Bunkerbestand": 0.5, + "Bunkerzufluss": 0.5, + "Vortags_Bunkerbestand": 0.5, + }, + index_scale=1.2, + ) + for i, w in enumerate(widths): + worksheet.set_column(i, i, w) + + header_fmt = workbook.add_format({"bold": True, "bg_color": "#E6E6E6", "border": 1}) + block_colors = ["#DCEFFE", "#FDEBD0", "#E8F8F5", "#FADBD8", "#E8DAEF", "#FEF9E7"] + block_formats = [ + workbook.add_format({"bold": True, "bg_color": color, "border": 1}) for color in block_colors + ] + + index_cols = len(bunker_sheet.index.names) + n_header_rows = bunker_sheet.columns.nlevels + + # Base header formatting for index columns. + for r in range(n_header_rows): + for c in range(index_cols): + worksheet.write(r, c, "", header_fmt) + + # Apply block colors per Kraftwerk on header rows. + top_level = bunker_sheet.columns.get_level_values(0) + blocks = {} + for idx, label in enumerate(top_level): + blocks.setdefault(label, []).append(idx) + + for b_idx, (label, cols) in enumerate(blocks.items()): + fmt = block_formats[b_idx % len(block_formats)] + for r in range(n_header_rows): + for c in cols: + value = bunker_sheet.columns.get_level_values(r)[c] + worksheet.write(r, index_cols + c, value, fmt) + else: + from openpyxl.styles import Border, Font, PatternFill, Side + + worksheet = writer.sheets["mit_Bunkerbestand"] + widths = autosize_worksheet(worksheet, bunker_sheet) + widths = adjust_widths_for_labels( + bunker_sheet, + widths, + { + "Reichwalde": 0.5, + "Nochten": 0.5, + "Welzow": 0.5, + "Gesamt": 0.5, + "Nachfrage": 0.5, + "Bunkerbestand": 0.5, + "Bunkerzufluss": 0.5, + "Vortags_Bunkerbestand": 0.5, + }, + index_scale=1.2, + ) + for i, w in enumerate(widths): + worksheet.column_dimensions[chr(65 + i)].width = w + header_fill = PatternFill("solid", fgColor="E6E6E6") + block_colors = ["DCEFFE", "FDEBD0", "E8F8F5", "FADBD8", "E8DAEF", "FEF9E7"] + block_fills = [PatternFill("solid", fgColor=c) for c in block_colors] + bold_font = Font(bold=True) + border = Border( + left=Side(style="thin"), + right=Side(style="thin"), + top=Side(style="thin"), + bottom=Side(style="thin"), + ) + + index_cols = len(bunker_sheet.index.names) + n_header_rows = bunker_sheet.columns.nlevels + top_level = bunker_sheet.columns.get_level_values(0) + blocks = {} + for idx, label in enumerate(top_level): + blocks.setdefault(label, []).append(idx) + + for r in range(n_header_rows): + for c in range(index_cols): + cell = worksheet.cell(row=r + 1, column=c + 1) + cell.fill = header_fill + cell.font = bold_font + cell.border = border + + for b_idx, (label, cols) in enumerate(blocks.items()): + fill = block_fills[b_idx % len(block_fills)] + for r in range(n_header_rows): + for c in cols: + cell = worksheet.cell(row=r + 1, column=index_cols + c + 1) + cell.fill = fill + cell.font = bold_font + cell.border = border + + # Kohlesorten-Mischverhaeltnis (gesamter Zeitraum) + j_name_map = { + "J": "Jänschwalde", + "SP": "Schwarze Pumpe", + "B3": "Boxberg Werk 3", + "B4": "Boxberg Werk 4", + } + i_name_map = { + "Reichwalde": "Reichwalder-Kohle", + "Nochten": "Nochtener-Kohle", + "Welzow": "Welzower-Kohle", + } + + # Empirical mix over full horizon based on delivered quantities (x). + total_delivered_by_j = {} + for j_code in j_name_map: + if j_code not in model.J: + continue + total_delivered_by_j[j_code] = sum( + safe_value(model.x[i, j_code, w, d, s]) + for i in model.I + for w in model.W + for d in model.D + for s in model.S + if (i, j_code, w, d, s) in model.x + ) + + # Empirical bunker mix over full horizon based on bunker stock. + total_bunker_by_j = {} + if hasattr(model, "bunker"): + for j_code in j_name_map: + if j_code not in getattr(model, "J_BUNKER", []): + continue + total_bunker_by_j[j_code] = sum( + safe_value(model.bunker[i, j_code, w, d]) + for i in model.I + for w in model.W + for d in model.D + if (i, j_code, w, d) in model.bunker + ) + + mix_rows = [] + for j_code, j_name in j_name_map.items(): + if j_code not in model.J: + continue + for i_code, i_name in i_name_map.items(): + if i_code not in model.I: + continue + denom = total_delivered_by_j.get(j_code, 0.0) + num = sum( + safe_value(model.x[i_code, j_code, w, d, s]) + for w in model.W + for d in model.D + for s in model.S + if (i_code, j_code, w, d, s) in model.x + ) + empirisch = round(100 * num / denom, 2) if denom > 0 else 0.0 + bunker_empirisch = 0.0 + if hasattr(model, "bunker") and j_code in total_bunker_by_j: + denom_b = total_bunker_by_j.get(j_code, 0.0) + num_b = sum( + safe_value(model.bunker[i_code, j_code, w, d]) + for w in model.W + for d in model.D + if (i_code, j_code, w, d) in model.bunker + ) + bunker_empirisch = round(100 * num_b / denom_b, 2) if denom_b > 0 else 0.0 + mix_rows.append( + { + "kraftwerk": j_name, + "kohlesorte": i_name, + "ziel_low": round(100 * pyo.value(model.alpha_target_low[i_code, j_code]), 2), + "ziel_high": round(100 * pyo.value(model.alpha_target_high[i_code, j_code]), 2), + "maximal": round(100 * pyo.value(model.alpha_max[i_code, j_code]), 2), + "minimal": round(100 * pyo.value(model.alpha_min[i_code, j_code]), 2), + "empirisch": empirisch, + "bunker_empirisch": bunker_empirisch, + } + ) + + mix_df = pd.DataFrame(mix_rows) + mix_df.to_excel(writer, sheet_name="Kohlemischverhältnis", index=False) + if excel_engine == "xlsxwriter": + ws_mix = writer.sheets["Kohlemischverhältnis"] + workbook = writer.book + mix_block_colors = ["#E8F8F5", "#FDEBD0", "#DCEFFE", "#FADBD8"] + mix_formats = [workbook.add_format({"bg_color": c}) for c in mix_block_colors] + red_fill = workbook.add_format({"bg_color": "#F5B7B1"}) + green_fill = workbook.add_format({"bg_color": "#D4EFDF"}) + if not mix_df.empty: + emp_col = mix_df.columns.get_loc("empirisch") + bunker_emp_col = mix_df.columns.get_loc("bunker_empirisch") + current = None + block_idx = -1 + for r, row in mix_df.iterrows(): + if row["kraftwerk"] != current: + block_idx += 1 + current = row["kraftwerk"] + fmt = mix_formats[block_idx % len(mix_formats)] + for c in range(0, min(6, mix_df.shape[1])): + ws_mix.write(r + 1, c, row.iloc[c], fmt) + emp_fmt = ( + red_fill + if row["empirisch"] < row["ziel_low"] or row["empirisch"] > row["ziel_high"] + else green_fill + ) + ws_mix.write(r + 1, emp_col, row["empirisch"], emp_fmt) + ws_mix.write(r + 1, bunker_emp_col, row["bunker_empirisch"]) + widths = autosize_worksheet(ws_mix, mix_df, index_cols=[]) + for i, w in enumerate(widths[1:]): + ws_mix.set_column(i, i, w) + else: + ws_mix = writer.sheets["Kohlemischverhältnis"] + if not mix_df.empty: + from openpyxl.styles import PatternFill + + mix_block_colors = ["E8F8F5", "FDEBD0", "DCEFFE", "FADBD8"] + mix_fills = [PatternFill("solid", fgColor=c) for c in mix_block_colors] + red_fill = PatternFill("solid", fgColor="F5B7B1") + green_fill = PatternFill("solid", fgColor="D4EFDF") + current = None + block_idx = -1 + for row_idx, kraftwerk in enumerate(mix_df["kraftwerk"], start=2): + if kraftwerk != current: + block_idx += 1 + current = kraftwerk + fill = mix_fills[block_idx % len(mix_fills)] + for col_idx in range(1, min(7, mix_df.shape[1]) + 1): + ws_mix.cell(row=row_idx, column=col_idx).fill = fill + emp_col = mix_df.columns.get_loc("empirisch") + 1 + bunker_emp_col = mix_df.columns.get_loc("bunker_empirisch") + 1 + for r, row in mix_df.iterrows(): + fill = red_fill if row["empirisch"] < row["ziel_low"] or row["empirisch"] > row["ziel_high"] else green_fill + ws_mix.cell(row=r + 2, column=emp_col).fill = fill + # No group fill for bunker_empirisch (column H). + widths = autosize_worksheet(ws_mix, mix_df, index_cols=[]) + for i, w in enumerate(widths[1:]): + ws_mix.column_dimensions[chr(65 + i)].width = w + + +def main() -> None: + parser = argparse.ArgumentParser(description="Run the Pyomo optimization model.") + parser.add_argument( + "--data-dir", + type=Path, + default=Path("data/processed"), + help="Directory containing input parquet files.", + ) + parser.add_argument( + "--solver", + default="gurobi", + help="Solver name passed to Pyomo (default: gurobi).", + ) + parser.add_argument( + "--max-rows", + type=int, + default=50, + help="Maximum number of non-zero decision variable rows to print.", + ) + parser.add_argument( + "--time-limit", + type=int, + default=600, + help="Time limit (seconds) for the solver (default: 600).", + ) + parser.add_argument( + "--output-xlsx", + type=Path, + default=Path("data/out/output.xlsx"), + help="Excel output file for deliveries by plant/week/day/shift.", + ) + parser.add_argument( + "--mip-gap", + type=float, + default=0.03, + help="MIP gap tolerance (default: 0.03).", + ) + parser.add_argument( + "--step-size-tonnes", + type=int, + default=1000, + choices=[960, 1000], + help="Discrete train step size in tonnes (default: 1000).", + ) + args = parser.parse_args() + + tables = load_tables(args.data_dir) + model = build_model(tables, step_size_tonnes=args.step_size_tonnes) + + solve_model(model, args.solver, args.time_limit, args.mip_gap) + # report_results(model, args.max_rows) + export_results(model, args.output_xlsx) + + +if __name__ == "__main__": + main() + + +# uv run python src/optimization/run_optimization.py --solver gurobi --mip-gap 0.05 +# uv run python src/optimization/run_optimization.py --solver highs diff --git a/src/preprocessing/exploration_preprocess.py b/src/preprocessing/exploration_preprocess.py new file mode 100644 index 0000000..3604515 --- /dev/null +++ b/src/preprocessing/exploration_preprocess.py @@ -0,0 +1,505 @@ +# Generated from exploration.ipynb +# %% +import os +from pathlib import Path + +import numpy as np +import pandas as pd + +# %% +PROJECT_ROOT = Path(__file__).resolve().parents[2] +DEFAULT_INPUT = PROJECT_ROOT / "data/input/PoC1_Rohkohleverteilung_Input_Parameter.xlsx" +DEFAULT_OUTPUT = PROJECT_ROOT / "data/processed" +INPUT_XLSX = Path(os.environ.get("POC1_INPUT_XLSX", str(DEFAULT_INPUT))) +OUTPUT_DIR = Path(os.environ.get("POC1_OUTPUT_DIR", str(DEFAULT_OUTPUT))) +OUTPUT_DIR.mkdir(parents=True, exist_ok=True) +path = INPUT_XLSX +# %% [markdown] +# # Mappe Parameter +# %% [markdown] +# ## Erlaubte Abweichung +# ### Kraftwerke +# %% +raw = pd.read_excel(path, sheet_name="Parameter", header=None) + +t1 = raw.iloc[0:18, 0:8].copy() +# 1. komplett leere Zeilen raus +df = t1.dropna(how="all").reset_index(drop=True) + +# 2. Einheitenzeile und Labelzeile identifizieren +unit_idx = df[df.apply(lambda r: r.astype(str).str.contains(r"\[kt\]").any(), axis=1)].index[0] +label_idx = unit_idx + 1 + +unit_row = df.loc[unit_idx] +label_row = df.loc[label_idx] + +# 3. Spaltennamen bauen +headers = [] +for u, l in zip(unit_row, label_row): + if pd.isna(l): + headers.append(str(u)) + elif pd.isna(u): + headers.append(str(l)) + else: + headers.append(f"{l} ({u})") + +# 4. Daten unterhalb der Headerzeilen +data = df.loc[label_idx+1:].reset_index(drop=True) +data.columns = headers + +# Neue Kopie mit eindeutigen Spaltennamen +data2 = data.copy() +data2.columns = [ + "col0", + "titel", + "col2", + "zeitraum", + "minus_kt", + "plus_kt", + "minus_pct", + "plus_pct", +] + +# Titelzeilen entfernen +mask_title = data2["titel"].isin([ + "Erlaubte Abweichungen der Bedarfserfüllung", + "Kraftwerk", + "Kraftwerke", + "Veredlung ISP" +]) +data3 = data2[~mask_title].reset_index(drop=True) + +# kraftwerk steht in col2, nach unten füllen +data3["kraftwerk"] = data3["col2"].ffill() + +# Zeilen ohne Zeitraum entfernen +data3 = data3.dropna(subset=["zeitraum"]).reset_index(drop=True) + +# numerische Spalten konvertieren +for col in ["minus_kt", "plus_kt", "minus_pct", "plus_pct"]: + data3[col] = pd.to_numeric(data3[col], errors="coerce") + +# Endauswahl +result = data3[["kraftwerk", "zeitraum", "minus_kt", "plus_kt", "minus_pct", "plus_pct"]] +data3 = data2.copy() + +data3["kraftwerk"] = data3["col2"].ffill() +data3 = data3.dropna(subset=["zeitraum"]).reset_index(drop=True) + +for col in ["minus_kt", "plus_kt", "minus_pct", "plus_pct"]: + data3[col] = pd.to_numeric(data3[col], errors="coerce") + +bounds_power_plants = data3[["kraftwerk", "zeitraum", "minus_kt", "plus_kt", "minus_pct", "plus_pct"]] +bounds_power_plants[["minus"]] = bounds_power_plants[["minus_kt"]]*1000 +bounds_power_plants[["plus"]] = bounds_power_plants[["plus_kt"]]*1000 +bounds_power_plants.drop(columns=["minus_kt", "plus_kt"],inplace=True) +bounds_power_plants.to_parquet(OUTPUT_DIR / "bounds_power_plants.parquet") +print("Saved bounds_power_plants.parquet") +bounds_power_plants +# %% [markdown] +# ### Veredlung +# %% +raw = pd.read_excel(path, sheet_name="Parameter", header=None) + +t1_upper = raw.iloc[0:4, 0:8].copy() +t1_upper +t2 = raw.iloc[18:28, 0:8].copy() + +t2 = pd.concat([t1_upper,t2], axis=0) +# 1. komplett leere Zeilen raus +df = t2.dropna(how="all").reset_index(drop=True) + +# 2. Einheitenzeile und Labelzeile identifizieren +unit_idx = df[df.apply(lambda r: r.astype(str).str.contains(r"\[kt\]").any(), axis=1)].index[0] +label_idx = unit_idx + 1 + +unit_row = df.loc[unit_idx] +label_row = df.loc[label_idx] + +# 3. Spaltennamen bauen +headers = [] +for u, l in zip(unit_row, label_row): + if pd.isna(l): + headers.append(str(u)) + elif pd.isna(u): + headers.append(str(l)) + else: + headers.append(f"{l} ({u})") + +# 4. Daten unterhalb der Headerzeilen +data = df.loc[label_idx+1:].reset_index(drop=True) +data.columns = headers + +# Neue Kopie mit eindeutigen Spaltennamen +data2 = data.copy() +data2.columns = [ + "col0", + "titel", + "col2", + "zeitraum", + "minus_kt", + "plus_kt", + "minus_pct", + "plus_pct", +] + +# Titelzeilen entfernen +mask_title = data2["titel"].isin([ + "Erlaubte Abweichungen der Bedarfserfüllung", + "Kraftwerk", + "Kraftwerke", + "Veredlung ISP" +]) +# kraftwerk steht in col2, nach unten füllen +data3["kraftwerk"] = data3["col2"].ffill() + +# Zeilen ohne Zeitraum entfernen +data3 = data3.dropna(subset=["zeitraum"]).reset_index(drop=True) + +# numerische Spalten konvertieren +for col in ["minus_kt", "plus_kt", "minus_pct", "plus_pct"]: + data3[col] = pd.to_numeric(data3[col], errors="coerce") + +# Endauswahl +result = data3[["kraftwerk", "zeitraum", "minus_kt", "plus_kt", "minus_pct", "plus_pct"]] +data3 = data2.copy() + +data3["kraftwerk"] = data3["col2"].ffill() +data3 = data3.dropna(subset=["zeitraum"]).reset_index(drop=True) + +for col in ["minus_kt", "plus_kt", "minus_pct", "plus_pct"]: + data3[col] = pd.to_numeric(data3[col], errors="coerce") + +veredelung_bounds = data3[["kraftwerk", "zeitraum", "minus_kt", "plus_kt", "minus_pct", "plus_pct"]] +# veredelung_bounds.rename({"kraftwerk":"kohleart"}, inplace=True) +veredelung_bounds = veredelung_bounds.rename(columns={"kraftwerk":"kohlesorte"}) +veredelung_bounds[["minus"]] = veredelung_bounds[["minus_kt"]]*1000 +veredelung_bounds[["plus"]] = veredelung_bounds[["plus_kt"]]*1000 +veredelung_bounds.drop(columns=["minus_kt", "plus_kt"],inplace=True) +veredelung_bounds.to_parquet(OUTPUT_DIR / "veredelung_bounds.parquet") +print("Saved veredelung_bounds.parquet") +veredelung_bounds +# %% [markdown] +# ## Kohlesorten-Mischverhältnis +# %% + + +raw = pd.read_excel(path, sheet_name="Parameter", header=None) + +# J3:P16 -> Zeilen 2:16, Spalten 9:16 (0-basiert, rechte Grenze exklusiv) +block = raw.iloc[2:16, 9:16].copy() + +# Leere Zeilen entfernen +df = block.dropna(how="all").reset_index(drop=True) + +# Header finden: Einheitenzeile enthält "[%]" +unit_idx = df[df[12].astype(str).str.contains(r"\[%\]", regex=True)].index[0] +label_idx = unit_idx + 1 + +unit_row = df.loc[unit_idx] +label_row = df.loc[label_idx] + +# Spaltenüberschriften konstruieren +headers = [] +for u, l in zip(unit_row, label_row): + if pd.isna(l): + headers.append(str(u)) + elif pd.isna(u): + headers.append(str(l)) + else: + headers.append(f"{l} ({u})") + +# Daten unterhalb der Headerzeilen +data = df.loc[label_idx+1:].reset_index(drop=True) +data.columns = headers + +# Spalten sinnvoll benennen +data2 = data.copy() +data2.columns = [ + "titel", + "kraftwerk", + "kohlesorte", + "ziel_low", + "ziel_high", + "maximal", + "minimal", +] + +# Kraftwerk nach unten füllen +data2["kraftwerk"] = data2["kraftwerk"].ffill() + +# *** WICHTIG: nach kohlesorte filtern, nicht nach titel *** +data3 = data2[data2["kohlesorte"].notna()].reset_index(drop=True) + +# Finaler DataFrame +kohle_mix = data3[[ + "kraftwerk", + "kohlesorte", + "ziel_low", + "ziel_high", + "maximal", + "minimal", +]] + +# numerische Spalten in konsistente floats umwandeln +num_cols = ["ziel_low", "ziel_high", "maximal", "minimal"] +kohle_mix[num_cols] = kohle_mix[num_cols].apply(pd.to_numeric, errors="coerce") + +kohle_mix.to_parquet(OUTPUT_DIR / "kohle_mix.parquet") +print("Saved kohle_mix.parquet") +kohle_mix + +# %% [markdown] +# ## Förderkapazitäten +# %% + +# J19:M23 -> Zeilen 18:23, Spalten 9:13 (0-basiert) +block = pd.read_excel(path, sheet_name="Parameter", header=None).iloc[18:23, 9:13].copy() + +# leere Zeilen entfernen +df = block.dropna(how="all").reset_index(drop=True) + +# Titelzeile "Förderkapazität" rausfiltern +df = df[df[9] != "Förderkapazität"].reset_index(drop=True) + +# Spalten benennen +df.columns = ["kategorie", "tagebau", "zeitraum", "maximal"] + +# "Tagebau" nach unten auffüllen (für Nochten, Gesamt, Welzow-Süd) +df["kategorie"] = df["kategorie"].ffill() + +# Maximalwert in Zahl konvertieren +df["maximal"] = pd.to_numeric(df["maximal"], errors="coerce") + +# falls du nur die wesentlichen Infos brauchst: +foerderkap = df[["tagebau", "zeitraum", "maximal"]] + +foerderkap["maximal"] = foerderkap["maximal"]*1000 +foerderkap.to_parquet(OUTPUT_DIR / "foerderkapaz.parquet") +print("Saved foerderkap.parquet") +foerderkap +# %% [markdown] +# ## Verladungskapazitäten +# %% +raw = pd.read_excel(path, sheet_name="Parameter", header=None) + +# J26:M30 -> rows 25:30, cols 9:13 +block = raw.iloc[25:30, 9:13].copy() + +# vollständig leere Zeilen entfernen +df = block.dropna(how="all").reset_index(drop=True) + +# Titelzeile entfernen +df = df[df[9] != "Verladungskapazität"].reset_index(drop=True) + +# Spalten sinnvoll benennen +df.columns = ["kategorie", "verladung", "zeitraum", "maximal"] + +# Verladung nach unten auffüllen +df["verladung"] = df["verladung"].ffill() + +# Maximalwert numerisch +df["maximal"] = pd.to_numeric(df["maximal"], errors="coerce") + +# finale Auswahl +verladung = df[["verladung", "zeitraum", "maximal"]] +verladung["maximal"] = verladung.maximal*1000 +verladung.to_parquet(OUTPUT_DIR / "verladungskap.parquet") +print("Saved verladungskap.parquet") +verladung +# %% [markdown] +# ## Zugdurchlass +# %% +raw = pd.read_excel(path, sheet_name="Parameter", header=None) +block = raw.iloc[3:21, 17:24].copy() +block.replace("unlimitiert", np.inf, inplace=True) +# komplett leere Zeilen entfernen +df = block.dropna(how="all").reset_index(drop=True) + +# Kopfzeilen / Überschrift entfernen +mask_header = df.apply( + lambda r: r.astype(str).str.contains( + "Zugdurchlasskapazität|Maximal|Vielfaches von", regex=True + ).any(), + axis=1, +) +df = df[~mask_header].reset_index(drop=True) + +# Spalten benennen +df.columns = ["von", "start", "zum", "ziel", "zeitraum", "maximal", "vielfaches_von"] + +# numerische Spalten nach float konvertieren +for col in ["maximal", "vielfaches_von"]: + df[col] = pd.to_numeric(df[col], errors="coerce") + +# alle numerischen Werte * 1000 +df["maximal"] = df["maximal"] * 1000 + +# fertiger DataFrame +zugdurchlass = df.copy() +zugdurchlass.to_parquet(OUTPUT_DIR / "zugdurchlass.parquet") +zugdurchlass +# %% [markdown] +# # Mappe Rohkohlebedarf +# %% +raw = pd.read_excel(path, sheet_name="Rohkohlebedarf", header=None) + +df = raw.iloc[2:36, 1:15].copy().reset_index(drop=True) + +jahr = int(df.loc[0, 2]) +monat = str(df.loc[1, 2]) + +kw_header = df.loc[1, 4:10].tolist() +ver_header = df.loc[1, 12:14].tolist() + +kw_names = kw_header[:-1] + ["Gesamt_KW"] +ver_names = ver_header[:-1] + ["Gesamt_Veredlung"] + +data = df.loc[3:].reset_index(drop=True) + +out = pd.DataFrame() +out["jahr"] = jahr +out["monat"] = monat +out["datum"] = pd.to_datetime(data[1]) + +for idx, name in zip(range(4, 4+len(kw_names)), kw_names): + out[name] = pd.to_numeric(data[idx], errors="coerce") + +for idx, name in zip(range(12, 12+len(ver_names)), ver_names): + out[name] = pd.to_numeric(data[idx], errors="coerce") + +# units: convert kt to t +kw_cols = kw_names +out[kw_cols] = out[kw_cols] * 1000 + +# convert Gesamt_Veredlung from kt to t +welz_col, nocht_col, ges_ver_col = ver_names +out[ges_ver_col] = out[ges_ver_col] * 1000 + +out = pd.DataFrame() +out["datum"] = pd.to_datetime(data[1]) +out["jahr"] = jahr +out["monat"] = monat + +for idx, name in zip(range(4, 4+len(kw_names)), kw_names): + out[name] = pd.to_numeric(data[idx], errors="coerce") + +for idx, name in zip(range(12, 12+len(ver_names)), ver_names): + out[name] = pd.to_numeric(data[idx], errors="coerce") + +kw_cols = kw_names +out[kw_cols] = out[kw_cols] * 1000 +welz_col, nocht_col, ges_ver_col = ver_names +out[ges_ver_col] = out[ges_ver_col] * 1000 + +out.rename(columns={"Welzower Kohle":"Veredel_Welzower", "Nochtener Kohle": "Veredel_Nochtener"}, inplace=True) +out.to_parquet(OUTPUT_DIR / "rohkohlebedarf.parquet") +print("Saved rohkohlebedarf.parquet") +out.round(5) + +# %% [markdown] +# # Mappe Verfügbarkeit +# %% +import pandas as pd + + +raw = pd.read_excel(path, sheet_name="Verfügbarkeit", header=None) + +# Jahr/Monat bleiben statisch in C3/C4 +jahr = int(raw.iloc[2, 2]) +monat = str(raw.iloc[3, 2]) + +# B8:J38 (Datum in Spalte B, Wochentag in Spalte C) +df = raw.iloc[7:38, 1:10].copy().reset_index(drop=True) + +# "Datum"-Zeile entfernen +data = df.copy() +data = data[data[1] != "Datum"].reset_index(drop=True) +data = data[data[1].notna()].reset_index(drop=True) + +ver = pd.DataFrame() +ver["datum"] = pd.to_datetime(data[1]) +ver["jahr"] = jahr +ver["monat"] = monat + +# Welzow-Süd Schicht 1 2 3 in cols 4,5,6 +ver["Welzow_Sued_S1_t"] = pd.to_numeric(data[4], errors="coerce") * 1000 +ver["Welzow_Sued_S2_t"] = pd.to_numeric(data[5], errors="coerce") * 1000 +ver["Welzow_Sued_S3_t"] = pd.to_numeric(data[6], errors="coerce") * 1000 + +# Boxberg (NO+RW) Schicht 1 2 3 in cols 7,8,9 +ver["Boxberg_NO_RW_S1_t"] = pd.to_numeric(data[7], errors="coerce") * 1000 +ver["Boxberg_NO_RW_S2_t"] = pd.to_numeric(data[8], errors="coerce") * 1000 +ver["Boxberg_NO_RW_S3_t"] = pd.to_numeric(data[9], errors="coerce") * 1000 + +# ver.to_parquet +ver.to_parquet(OUTPUT_DIR / "Verfuegbarkeiten.parquet") +print("Saved Verfuegbarkeiten.parquet") +ver.round(5) + +# KVB Nord Zugdurchlasskapazitäten (L8:N38) +kvb_block = raw.iloc[7:38, 11:14].copy().reset_index(drop=True) +kvb_block = kvb_block.iloc[: len(data)].reset_index(drop=True) + +kvb = pd.DataFrame() +kvb["datum"] = pd.to_datetime(data[1]) +kvb["jahr"] = jahr +kvb["monat"] = monat +kvb["KVB_Nord_S1_t"] = pd.to_numeric(kvb_block[11], errors="coerce") * 1000 +kvb["KVB_Nord_S2_t"] = pd.to_numeric(kvb_block[12], errors="coerce") * 1000 +kvb["KVB_Nord_S3_t"] = pd.to_numeric(kvb_block[13], errors="coerce") * 1000 + +kvb.to_parquet(OUTPUT_DIR / "zugdurchlass_kvb_nord.parquet") +print("Saved zugdurchlass_kvb_nord.parquet") + + +# %% [markdown] +# # Bunker + +raw = pd.read_excel(path, sheet_name="Parameter", header=None) + +# Ausschnitt R23:W30 (0-basiert: Zeilen 22–29, Spalten 17–22) +bunker = raw.iloc[22:30, 17:23].copy() + +# Kraftwerks- und Veredlungsbunker (Jänschwalde, SP, BW3, ISP) +plants = bunker.iloc[25-22:29-22].reset_index(drop=True) +plants.columns = [ + "typ", + "anlage", + "anfang_mo_di_kt", + "anfang_rest_kt", + "zielbestand_kt", + "maximal_kt", +] + +# Typ auffüllen +plants["typ"] = plants["typ"].ffill() + +# numerische Werte konvertieren +for col in ["anfang_mo_di_kt", "anfang_rest_kt", "zielbestand_kt", "maximal_kt"]: + plants[col] = pd.to_numeric(plants[col], errors="coerce") + +# *** Umrechnung in t *** +plants = plants.rename(columns={ + "anfang_mo_di_kt": "anfang_mo_di_t", + "anfang_rest_kt": "anfang_rest_t", + "zielbestand_kt": "zielbestand_t", + "maximal_kt": "maximal_t", +}) + +plants[["anfang_mo_di_t", "anfang_rest_t", "zielbestand_t", "maximal_t"]] *= 1000 + +# Vorfahrfenster +vorfahrfenster_tage = pd.to_numeric(bunker.loc[29, 19], errors="coerce") + +plants.to_parquet(OUTPUT_DIR / "bunker.parquet") +print("Saved bunker.parquet") +plants + +pd.DataFrame([{"vorfahrfenster_tage": vorfahrfenster_tage}]).to_parquet( + OUTPUT_DIR / "bunker_vorfahrfenster.parquet" +) +print("Saved bunker_vorfahrfenster.parquet") + +# %% +print("\n ####################### Done preprocessing exploration_preprocess.py ####################### \n")