From cfe1e3035746454adfc02f66c6190fe58bb7101e Mon Sep 17 00:00:00 2001 From: Ryohei Sasaki Date: Mon, 22 Jan 2024 22:02:22 +0900 Subject: [PATCH] Addition of velocity scale factor estimation to EKF (#937) * Addition of velocity scale factor estimation to EKF * Format * Add a scale factor motion model in the Jacobian function description * Fix Jacobian function description * New script: 'ekf_with_velocity_correction.py' * Add doc * Fix doc * Correct the parts where the first letter of the sentence is lowercase * Fix doc title * Fix script title * Do grouping * Fix wrong motion function in ekf doc * Update docs/modules/localization/extended_kalman_filter_localization_files/extended_kalman_filter_localization_main.rst --------- Co-authored-by: Atsushi Sakai --- .../ekf_with_velocity_correction.py | 198 ++++++++++++++++++ .../ekf_with_velocity_correction_1_0.png | Bin 0 -> 40148 bytes ...tended_kalman_filter_localization_main.rst | 105 ++++++++++ 3 files changed, 303 insertions(+) create mode 100644 Localization/extended_kalman_filter/ekf_with_velocity_correction.py create mode 100644 docs/modules/localization/extended_kalman_filter_localization_files/ekf_with_velocity_correction_1_0.png diff --git a/Localization/extended_kalman_filter/ekf_with_velocity_correction.py b/Localization/extended_kalman_filter/ekf_with_velocity_correction.py new file mode 100644 index 00000000..5dd97830 --- /dev/null +++ b/Localization/extended_kalman_filter/ekf_with_velocity_correction.py @@ -0,0 +1,198 @@ +""" + +Extended kalman filter (EKF) localization with velocity correction sample + +author: Atsushi Sakai (@Atsushi_twi) +modified by: Ryohei Sasaki (@rsasaki0109) + +""" +import sys +import pathlib + +sys.path.append(str(pathlib.Path(__file__).parent.parent.parent)) + +import math +import matplotlib.pyplot as plt +import numpy as np + +from utils.plot import plot_covariance_ellipse + +# Covariance for EKF simulation +Q = np.diag([ + 0.1, # variance of location on x-axis + 0.1, # variance of location on y-axis + np.deg2rad(1.0), # variance of yaw angle + 0.4, # variance of velocity + 0.1 # variance of scale factor +]) ** 2 # predict state covariance +R = np.diag([0.1, 0.1]) ** 2 # Observation x,y position covariance + +# Simulation parameter +INPUT_NOISE = np.diag([0.1, np.deg2rad(5.0)]) ** 2 +GPS_NOISE = np.diag([0.05, 0.05]) ** 2 + +DT = 0.1 # time tick [s] +SIM_TIME = 50.0 # simulation time [s] + +show_animation = True + + +def calc_input(): + v = 1.0 # [m/s] + yawrate = 0.1 # [rad/s] + u = np.array([[v], [yawrate]]) + return u + + +def observation(xTrue, xd, u): + xTrue = motion_model(xTrue, u) + + # add noise to gps x-y + z = observation_model(xTrue) + GPS_NOISE @ np.random.randn(2, 1) + + # add noise to input + ud = u + INPUT_NOISE @ np.random.randn(2, 1) + + xd = motion_model(xd, ud) + + return xTrue, z, xd, ud + + +def motion_model(x, u): + F = np.array([[1.0, 0, 0, 0, 0], + [0, 1.0, 0, 0, 0], + [0, 0, 1.0, 0, 0], + [0, 0, 0, 0, 0], + [0, 0, 0, 0, 1.0]]) + + B = np.array([[DT * math.cos(x[2, 0]) * x[4, 0], 0], + [DT * math.sin(x[2, 0]) * x[4, 0], 0], + [0.0, DT], + [1.0, 0.0], + [0.0, 0.0]]) + + x = F @ x + B @ u + + return x + + +def observation_model(x): + H = np.array([ + [1, 0, 0, 0, 0], + [0, 1, 0, 0, 0] + ]) + z = H @ x + + return z + + +def jacob_f(x, u): + """ + Jacobian of Motion Model + + motion model + x_{t+1} = x_t+v*s*dt*cos(yaw) + y_{t+1} = y_t+v*s*dt*sin(yaw) + yaw_{t+1} = yaw_t+omega*dt + v_{t+1} = v{t} + s_{t+1} = s{t} + so + dx/dyaw = -v*s*dt*sin(yaw) + dx/dv = dt*s*cos(yaw) + dx/ds = dt*v*cos(yaw) + dy/dyaw = v*s*dt*cos(yaw) + dy/dv = dt*s*sin(yaw) + dy/ds = dt*v*sin(yaw) + """ + yaw = x[2, 0] + v = u[0, 0] + s = x[4, 0] + jF = np.array([ + [1.0, 0.0, -DT * v * s * math.sin(yaw), DT * s * math.cos(yaw), DT * v * math.cos(yaw)], + [0.0, 1.0, DT * v * s * math.cos(yaw), DT * s * math.sin(yaw), DT * v * math.sin(yaw)], + [0.0, 0.0, 1.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 1.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 1.0]]) + return jF + + +def jacob_h(): + jH = np.array([[1, 0, 0, 0, 0], + [0, 1, 0, 0, 0]]) + return jH + + +def ekf_estimation(xEst, PEst, z, u): + # Predict + xPred = motion_model(xEst, u) + jF = jacob_f(xEst, u) + PPred = jF @ PEst @ jF.T + Q + + # Update + jH = jacob_h() + zPred = observation_model(xPred) + y = z - zPred + S = jH @ PPred @ jH.T + R + K = PPred @ jH.T @ np.linalg.inv(S) + xEst = xPred + K @ y + PEst = (np.eye(len(xEst)) - K @ jH) @ PPred + return xEst, PEst + + +def main(): + print(__file__ + " start!!") + + time = 0.0 + + # State Vector [x y yaw v s]' + xEst = np.zeros((5, 1)) + xEst[4, 0] = 1.0 # Initial scale factor + xTrue = np.zeros((5, 1)) + true_scale_factor = 0.9 # True scale factor + xTrue[4, 0] = true_scale_factor + PEst = np.eye(5) + + xDR = np.zeros((5, 1)) # Dead reckoning + + # history + hxEst = xEst + hxTrue = xTrue + hxDR = xTrue + hz = np.zeros((2, 1)) + + while SIM_TIME >= time: + time += DT + u = calc_input() + + xTrue, z, xDR, ud = observation(xTrue, xDR, u) + + xEst, PEst = ekf_estimation(xEst, PEst, z, ud) + + # store data history + hxEst = np.hstack((hxEst, xEst)) + hxDR = np.hstack((hxDR, xDR)) + hxTrue = np.hstack((hxTrue, xTrue)) + hz = np.hstack((hz, z)) + estimated_scale_factor = hxEst[4, -1] + if show_animation: + plt.cla() + # for stopping simulation with the esc key. + plt.gcf().canvas.mpl_connect('key_release_event', + lambda event: [exit(0) if event.key == 'escape' else None]) + plt.plot(hz[0, :], hz[1, :], ".g") + plt.plot(hxTrue[0, :].flatten(), + hxTrue[1, :].flatten(), "-b") + plt.plot(hxDR[0, :].flatten(), + hxDR[1, :].flatten(), "-k") + plt.plot(hxEst[0, :].flatten(), + hxEst[1, :].flatten(), "-r") + plt.text(0.45, 0.85, f"True Velocity Scale Factor: {true_scale_factor:.2f}", ha='left', va='top', transform=plt.gca().transAxes) + plt.text(0.45, 0.95, f"Estimated Velocity Scale Factor: {estimated_scale_factor:.2f}", ha='left', va='top', transform=plt.gca().transAxes) + plot_covariance_ellipse(xEst[0, 0], xEst[1, 0], PEst) + plt.axis("equal") + plt.grid(True) + plt.pause(0.001) + + +if __name__ == '__main__': + main() diff --git a/docs/modules/localization/extended_kalman_filter_localization_files/ekf_with_velocity_correction_1_0.png b/docs/modules/localization/extended_kalman_filter_localization_files/ekf_with_velocity_correction_1_0.png new file mode 100644 index 0000000000000000000000000000000000000000..595b651bd2354f67fd57de47e35b9848433d7c93 GIT binary patch literal 40148 zcmeFZWmr{R^gW8GbSWqeBGL%bf`lR=-5@O;(%qmaDk%ceC0!!jA>BwzcXu2byz{)j z|L;D}{eC~(5BGyF?{heN@3q&OYtAvp7>h58@{+h%WLPLDD7aFu#FSA`Q2kI)Zp>kz z!6!TeQ%mrJ-%0$nlZvgWldFM)35uM7lbw~Vla+_~$FF`l2-TIR&pZesBz2%)=vyDh)ESe7lq8J8zN9BEf@IN#k zI&|-%{P*EsQcRSO$d9&yZ}}mAL*OET`uxAIKcl^o^52&e!%-BGziH}wg~IsXjm#dQ z5Fx+(i|hZt&Hu;ehQd1W-DQszujBQlyP%NJz`9d;Unp`}sopCGZjKc$E-w1(*(N(1 z+(j|{eRjCIe|T8p%!^#*aLj>$fz$Kz2?26mnj3HKyLot^PT(#rFPkh*E$Mc4cIN9h zB*`ZV7(5}dN~ElLU}8hAyi=vnC3tf011Dds?CIH=W<&dS_$@!pHx|c;gM)7OZ1^m)FSW1ArcSEQexc)}hY^1AU$C3M(`gyO5KIvK ztPt4V-09!(_tjeh*HxqLO_nB&Cy$*bge<)6HvB!0*PT7kM##l#4z>r=!}Y4ad-rzn zMx-cVgwrDT%HNOAZr0ZL|4bQr;Z?rPMTCX&NtMyUp@roqLDE>=-sh7>uHU%*hb71%5IuF}?LH^oHxoa_%KRih zhAMcw)o41~qb)v{cxwN8?Bl{dnvcBs=hPjiK`3A z;^VZP6?8IXq~4b69m;THZ40?QMy7-~ch03x$|5OhT))VSe}9v{_VZYkz3y(8GrHLV zB}LWFnpmy>yG=C71xiXOyI}*DSK?VR;^=FU%aX1aRd~^=TzXQbSab|3jC7cVkz!aG zGwEN?(DYWw-R|J~p<@bu_lRp(ey9DM^^I(Np?MugWFv2y~Y!+p2~c3v8#B~gy$Dsyp=xw4@U??|-=hcYIz zx8s_8c0eq}u3E1=SyWar4LUphutr260^D=&)p@!3wl@aXFo#5s0q8~&W>Y5L$mK2 zkJip9hkEZ2RZyVj`mZ zSv9bxXwtV`?I-=`d|RAA?$!Y7ZGoGceqsqvH|F~9RD7dg+3G1$?xZ-Wllb7v&4(@C z++%Kz%3ZhkzVo@HB2^K7ln8pFKW}Ao1Eq;z*Dw1`iHS0+_Nr4ylt^|RVQa$D*6Ybw z3T#u!&n#pzIs>mCGMjl7i2n0X5-KL*3I3Sv{yUtkh%eio3vJlsh>@zH%)k9H$FA_J zw(#d!4o9y9*N(0f1)hNw-vdMt>_fcSdHeBa&T}EF{@(uwFVJ$ z8Ma^;8X1)&9Ct;roGk>gKXut6F)fKLHt8PBkRi@gNKRTVE$Lop4RT-4h~60&@xAua z^G1xYW(XNH-BK^qlY9O8_4h${qAy>*1XJ*TOZGYc63=HlSnut5ebVr6elvbwDC*m{ zrVp5eudAJyh}m>foYn_xPZxrMsYR&j^x`vi;NHQcTmgv9F+y^3nLPD^(~Z2sY@IV< zU&e)0WNzqi+`&0t*^Za&WYKe)(EUO!Ka$bfz!9ksaszMg;T$e2s?9(Cr|ee*W0*q-cw^p|7R^rlU^e-z_@!=$kP^@X9ia4zi-WO4iBmAM3g)RGS(0|C&ES zTiJ~dBQ4=ctb9l>8We+35Y0QgLlZhRVZKsEhMR5w`%BzEdPOd#zy+NQ?Lj6Flrh>q zJ37?kXX

VivrrRlmPaNRCsNQcs9hX)fpUe5e}G*;-oqP*LGpo<&YX)H0+X(wmj+ zBi)~btRxit4j*vId5tg1io3eIgBII|X5zp5J!O9-FW>g&iQsS%twiu#n8@|qYO-(U z>mtM2i=7&$-I>HX&y$6@xvx9?Dw(ndJ5yC&|E^Akvz7dXJdgJVW!ROn6=?$SsY8_m z8Is&k{eL&^Hv9{Xjwa?Z?MZq02@jrn%hmN{%d6&Ca``5#s`_+|t44-Qw7Q>GnWb9X zSe4M;&Q3CW{prJIZsNyUkvsf+wliL9UahUIOwZvzCdUQiAJfwxm3bg6qaq`XH}n67 zT&k}O*bi7Pu7_Ag+eed~^}jM9`D8g#mm8;39y{?i`v+0_nUr`|6Cz60G3#4yt(lVA zNvCwzzBhX7t(_~2kS6SpH@oVxYc~EV`<5;2VX4)%`5_*pnkxD+-~D1O zEll{zb-NoA&^6PJtS6pvz7la69k5|xYNW|FH?NEX%|$=6h= zs;Xkr*Vji(ye-Uua?O$vIZ)@B{8tlPGaQp+6sdv9{*gA&l2IVlu71-y%vEtY$y-)PJF9?fK~rV@x%iA)D7L*K znvjqvq9)ztJD85u?WexoSz@3I6D z(v>n%CE@Eg%^d^VFTN$%+z|W?T0_LYP|8n_7%PqTV+S~>}a)fLW(gdU->EM@W$RRNK6Bk!3L4PgThB@!%va&^IlCgRPrec|W#DSo##ueIPJUY! z5So)v;1F4&7W$qv<_1q*^`Hv%^NS%5flJQ->YmkgTv?%tkM^r`SRNi#8mo4=X76RF zOjg!ndjz>7XcX{LhnJNtx{P;c>W+42lQ9YD_l8vz?%ciG^vj@`Q9YjonLSY0Yxhu7 zQ&YXpmZNXox7vBNYQ+(paF`~3?cL$T+H4H>T$Mk7qvc!RG6+HLdApK>%yOsm^3VaFV;4}G#AfTQ4tCAsFX53dP zEk|<2MMY7meNJg0vC^}$;?(IWczFrJhRl?W<;>SCL4yGuZ%=+2w_3@`N?O>PZ_d}Q z2u(?O@H|gVwc_p2)b#XttveSJ6B7O7_7syDZbU5)UYqU)QOn6=AvLkV76ASD&y6sAQgLkdw9NR&@} z=E~N44O5y|WYo@_5$UlKYo;2_raM@B4Pm4?6c`$$9{-RpN4HJaQe}xA;q&K0TI5Bc zdFT&F}-#R@F7;S)VX6n6I4XI&z9{5KX9b<~cTQ@ik(Q%SNZ|Fw1T4YMt> zjntcYjucyWbDxp$SVrdR)m1BRG)%-dHa0p*U(aZj1`9{PtmW!?Zu4v&-l3I9&zY~e z&LIbvhZj#G0BKM3!of@;a3#1nBK-D8Eb@k#T&sNe4$qsv3Z^PTXYV zW)epW^$RH|C{$HdBc1i#+$vFrf9lrcuae@RIj;|9aIe;%%y?Qh`rQ~^b&Wm)1P8Ta zNa$?sZ_%RAvI9fOc>av%u29C%dSqH!+Tgm=2H*@i$3%zat{k7M^PIP1KXY1I#1P|V z$xP3kH-?#-ou~u6d{-mL8YS)b7sRj6H`(j|U0ial9S7T&J?+K;HCM9 z9hg2`ReBOg?#qa;QR^gHXHlhTF5G#Oh4?!i%5yrM#qEo(W8m9J6S?@E+hYCCzDhuUBtzzFWpL8!$#zFJYJW2CbC5 zv9xQba>YADLFAgK_RaPo!)x-JJ)iF5c~JDuWfKW}laj6z)g|0Yz@O!hG$r83zcj+s zr#>fhdG|9e)PiP+WZdczN^d-$*=W9IdS>R``}YI<{a?mu4fK$3n-j1%Tyg^5YW)5C z9U!Y2*C(H9YM%LCA&9f&6AfWhTYroEU1zX{eH0yMaFCIH#+hWEE+%4t9*LncrdZC(_qd4^r{!?=i`TD7 z0GKwmw|{zf6W#C$2_!hW!rB8YD9EkwzIc=Ze_-gc%aZA8bEru0Nu{8maDKGzIJ*d$ zx(ynG9JRdru#WKZL$&UXcoa{+&aeG|2z?vH?bm6HS{hmm^vL)LJr;?~_$$)oq*}!l ziHPt*SKr>?M6M@~$kyts*nBCnC09Q@p%_OqA1boqyY;Ks%U+#_Urx2v?t!43fPho2 z)~edXc}0d)EdA(QrhKQs%IuL;_$HQk;j$LiH(B;C$tmo z{rzG5jw_1QP8L0;TmWvC+d{};5*D|%f|Z$@i{<3x>Q)i3kDUajg;p6e@=e0y;$Rm{ z-W1{6VLOBniGU&xJBLv@{WGF6h?sq_(tZ(Vh@8igC`}@y*kvqNy?`uBA$bKd-7@60 z+bJ&F6I(NN@s*WNzkU10sPy9&?9D7lHk#*N5xZ%qBq(e;m0wqrJa^uCpFw33ji8sy z?}zS3zRm-IY)wxOy4i<>R!jPYjTTzIoNVbg_z2n9*l2cFQexoXlsezvlDCOIi%Mr{ zdDLb%cYedC?vzEX`R#sd`7utdo!_^uVycR6#YB^9I?b5_BI3)Y9VyAMq)_o^%H1XR zN2c2)8`Hj#m@R-Vz=lcn0h(Zbx$@%^2Ef&Zvc#W2;Z|9%EX02vJ9V&UaLTh)85?d zzqe<1mymEM{SGm)q@^YE+tK_~i4d~&`e@OQSo_dTDi_Z9zI+FlCDQkkLQg>jNU3#f zEE;>2{ZFyFjPx{XtH$F7!=05h7;)(g=` zU`BKr*UbR(dFp$8Ddf(smP@2vZj}nNcarW3Yvy$5bmDDlz!{4qJ8q)k%t0?9CL!Ul z)N$u|o&`4vPE4=qF{IFg&(ol9R{NZ^N89Lel~j&~x)5G&Gc+jV%;PW)Vb` z(=OtoDAFC^A$%^IC=hl3u@;Km_5@E;bMyL8mO{Wf^uzW`9pT)C$%5|R0Ey%D^ZQ;4 zdY_3{Sw5J_fy-#$T8D*;4nWa69+U)FLIn3*)USSqMd@(UGW4E5aMMmzc ziLUwma`c2gK0Zzx!03}A5L4wxrI2#0*UK_<8X71)u~le&c#qnW>>Era0Obts z-voYXX|%u-twKlffQKgvlHS6~N?R!9GdpKz2>eP)xRBq5)JcMa zR_(lr)|e#Z31v1LQMFuMR8(Yf+#AE;52*lkuZ2)y53nv0?E_YqfAtCjHW^w9fL$Us zotTH>flnM)WOa0OvUO{CVSiOtRTWuJ@~-qJ_eV0RWUFN1K~$PRg+!ua071&o%>pga z8cr*@HeSN;;>8OJ0VhVeB(A|sIr97LdOc7@B>{&G_PjI`v>g69`39?*np%D>NAY^< zlz}E%v{JpjW1?K*lnOP8sk!k_;E{;ud*KaK>z_MscC=0xcSUR4omI74OG}3Sq%IcW zg>jrrR(bJn%2Z_WTWJFAcD534t0=ucs%@(TgbH*A>r+)+0PItI|6Ql2r|3 zH9)bO=y;_o^N{K?#J=1GiNpSM&j8viQDDSK`D~`VccDg?&W-cLBZPY%?QmfjL!HWn zG#)EGX2qL_hD(+mt>>|x+8-{Tj?Ykt_lmK@A+HBX6CE^nk2ta*5Uhis0Qv+d6nL^RFxkK8Lnw+h%$Q!Hz-GycrxjPXIN7hv!3vdavArNO#mQ$5Q z2L}hrKrZFNR53F$c6j)hn`h_h)|5NAxTt>*YYxtQS9MVhFDYG82&HYKmv6T8`vs7l zs<5}B^78UFpc6Wu|3tkLD_{%+RW9r$i_CXY($dP%WTql}Bu)~XuE|Ml#Ll!frJ&m% zj~Qq${$qZeoz(yiD?i|l+?Hzo-hpK! zlL9w1*X!efB_KcPoTDsdsX{?kc!~w(8;7i#T|@{uUoJzkm-5-!8508ogIcaC(+H5p zW$>=ts|7k$$|E@{lKlMqIZ&qaguTyS>q8ojMZOktg5=_P)u7%e`AszV`Ubeqmc5PM z_*;~QtAB^QNLbIMM&j12p^j{_9YCtHa zfNXIn1&)%q_1xxfcSSN8ZTAS+FNh&A2BZs)SXm_hae2BJW(8HO*Ejmah%=T^HJiZR z-d_2p8})b7PypsqxNq1x*E#v(f_r8-;@2lLhpdfMBzC`KsV(@*m1G_Hux7sw*nJk@ zSij`FdHoFCeoE+z!wbbfyUuBy9i49(Fa__O&N^}V9Xhk0qmP)mx#1vGLN&wiU4Ax4 zBXnL5*%0TTl^kE#p0hQ_!U4FhhnG%64?cO3@C-XL?BZyv{U;oO;tIrvnA4h=t%HVx*4O<=}Um0_29B1 zX4cp&R{QORXx8GL&DInV^1b66+EZq0Rdb4@9t7?im$8VJ$f~a^I-o2>;;B$=1l2$I z)$TcM&Uu%NvMXSSJpV1v%}9QYz_C#^lUGLufoXEX-XTq1y1ahuo>X&|`>&Jo5EMdr zN9mg5rli+KDMRawdbLkW%=*hiaq1$3zE>51t{c6(>!o~?@x!l>oXIRsdUKo%HX?(A zINn+pAL@?-qpCuMES>pXQ88oQBJ_pUqmE8%ZVk&9PT7A_8+rBvJFm18%{k`<&r!w7&e``nh3F#sWA39dH$VB%VRdvASaFB zq1hU>*JpyJ84z&;uK6Q#M7+Ct+~xSW_>ziN`IKUwh6I@GnDHn>LR^as2-sKzvZ zH``%K`^Do~pwu?~K&BTt=Ji%E9ZABZ5u4epaCb*0t(23Dslo2{u-iWz1zzhWy>`n= zNB;_=atDv=2CoZeN3RgNnEL`D#a!K-kB#FBW&;RZURZ+irVbsQ0RtUMQWv+TP)Khd>yFMz%;c=?v=dExsp;s2D?<^zr!Bcw zS#V?5Rv}f+ki`F|m06=}-(`Q+?c`(rmh48QmYF~kW$La`gO-(vTm$3R3&-%l4te6dq%+nrV~VwdUKl zcT}B8F|G0qiy6a@5O0}eESspoaUx?F$ z;QuM7ahe$X+F!1lrPswF{AH?Zb4X^db9SD4cjSXHL&;?cjTqina7;Wi;UQ zk{6**tZV2wXRe8XdKo{fT(hV#7wN5sA1LFhy_pMUEv4JdPD-q|4agX%X{cP=G|Z(k zl&*G`A4%YVajEn1$Bq-w)0^ZW2-QCn&YflS6J9hWfAIBQG>qWS5&;s7eDn)zOqHt6_;cdxw?mva+Z8w*af`aQ3-YSG4-rhe3YjB zcKSDqtHLhRsObgXwc8hq!}ubrLwruH*PpE($Tf-**?qr3H6)rf;G}^vdmZt-SJm30 z^?9^M?ew9*{78jXShP}I)BT+->Zi0$XqX?3&dm2jSVOaWg4O8U79=osld$JwN$K3{ z#Caq*0wXjOI5Wpdw9Rxsdr5)80_)fR0pm$ySckoXZ`p|KqNrQa^S7<-J2vb?ds8u3r+h`%{HIlYC));?Gnj zRmuMuUL!+yiG1)DtJsaQrOUfBlY3NkW93V0(egMo)tEwow_D>6{%Uu!ZeHu}F9?y-zeuX9Pw79T^ ztw}RdEnmhjHhgDq{-BWwYJQmx9bpRovWsVBwn*%p?KGjFPIzbG4oJLWH#nUyR4+{< z>2`+fnpQIp#4g#sJ=wEATe%e+YU-Gl_I$!*GJRBG$r<>$o^|1N0mY>^1xDPYj=>dq zr8GDACT6e}=kTLRt7>VgRa;?KBUcq2Ho#}p#$uP$K(`fHdwVc@o z3hKJbL@_h2{xOfFd)^ulsYEKz*{%G( zcv{pS2BMDwsz|u)v?h7>%XU&<+C<6i<*1Fs&+UojeE5!glSJn~ZT&>^a48z>YYR#q zizOjlVP}6&Ll#N@tzLKuWnuc-%t%_W*|(&k%AKnkH_2#Axk#5CA4c3G5_`d!d&c-< zCe(18u-jZW{~q>5Vb%_p@RmdVi#rRabjJAD4ORi3wL!8E9*YUUrg}(S%zegow(nUk zSJ?9_C^cq)qqpWULV|2u`{EXg{|3%CADmX7}c3o_Pepy>f%E#<-O`=oc5aS1=ZzFci=C`yd6q=XPbW$0owfIbj|HwV6rcZs@ z-Z*&kEOU}?aEfKC>6zT@JcNgP**J^W-Uj`I%g}a_a3@xo}-7p zk#+xAgx<>R*DzJ{EA-8?G{^Z{~G+h9`9cB0*H$rX0~#ooLTvpw`d zoi^-_E-%_2G>?u1pQahj(fijGOFSOIQ?1!pl&{458{ejW=wp;@RWZqsdi_@jc0(y1c{J)T^!=q9uoJ9euBU`;scl=$yZW(`nMkXDsQ zxj7zs4^!rw#j99j*(dUlj$f-m4^zLFM9MhHXXLJ!kqX6oizdgxmqs zt1goQd7?{X$gvr*%GK*%am0(54V(O3E{GY}x5+5Sp{V~zo`mUl|57T#_P0{KmphO& zQD#Hka{S@0aPz!g{gpC4-hKI~Q3$G|uKU+>M(c!(TLEBMTB5-|ew;Fu z?{Hb?!86&wEFK;2Db*Mv{CI2OE7HOqM{6%Yq#N{p%*{=d zDBoW1ec?dI?mAKD$p?K;x$xM!Q@&oEK=*o*uy-6FPc+czfj5AbH3KLD5Xc^OTbN}X z?R?3q5h~!pte|S|G2=TF+pKO`l+V)k#5mJ8=t^e0XcT^D)|c=iYv|7X``iDxPbvY8YGXLZA82@h>_!{Ip7B-Sd60!(8c;fz53qU?pYbsivmZQ= z)HZvLk;)MljB{>8;EbofP^e}fsh;&T>oBvOQ2a0TFZ4g>7a_JhdU0OwVq!=(7eUQE zgkI1A#OqCA4UIAGcIc^$o{$s^Q6oun@AISit??3Ij?&$&p|Jt-{*wn{!{ox5IvqOX z@g!4a<-3&6T+*RAy*!?@jbzg$OB?8zD0}-QD5(2*(*)^{SR2d;LR3XW5T@{85<9T9 zkhJ|7oD36DcT=493j9MIIE47JA!}NUs7yHebz*F8&Zr=SAWHT+zyy`VbR<^|so_Fj zGXNZ`Q~w!o67D-yt4LB5$pTA;Q7w;X=(moHMDfkI-Sh+T8W^Cqlo$T5i;XcSD{Rf6 z*9YC;J@ky6+Ov5b{*fJ}K*Qe7hW}`r(J9D4+3_dtQ%=z z!pguWTJ$CIJ7$7FrxWLmP0XeW@1X);$%Pf&h_|ip+?xOq0rW~`ORyBTpE~?Jz))=-0&>pSFc%;x0OYtP%B}N&TZvRAo!vwlct(WWYaEsXJ<>HLxPWqC8F9m-L=(lP} z!-+{(q!Mtl%nu$sN7_L~H%*-sgoCfIF7nuPtJBBFkn{s+a%OIQT>XM+=f~Uern4Y$ zM}dI_;ra}m3&2%?V*-biQ;f%AaH;-`2_u>P;Hn^!`v!J6dhqO1szXUnO+tzQbj(}T z!)eR$#2n{|(e$&dBwVIAwsna65O^8DOXEBDU7j5RIYbe`ps=zwkmd!ih7`Dvwx@Nf z-L@G-S!?p7@19$TaU)Y)G%eZ zEqr09blD~YB9H#bldr%|?jx1wu(Qp=hP9zAd~jVog~c*!4#4Gx$(}xP*&RjWW7R7C z44c|zcSa64-1R@HH$lIVcHDw3=(gWV43@5wbQM?^Fmg0ppR?rT%;J{?et=wPOkB(> zvN*A@A|Xbg$x(nHMSZ}?h#b7n#Wo*`X@62LI4Y3LJUTvAI?VFX)_BC{&!4%+8ZMG> z@bDzEhU`~*VsQ>B$jGEYvmD$4d7q1hhK4jLA|~c@P>@WV_s^d{NkSz0V2{i?ootQ2 zbZ|KEmR2Id;Kh4>$a99jFi-U1-+f1EC!JEVgLVR*OiS8~=S+&&L8YppMwIsm6|a2P z3x6#Bpyt>0WH}7BL@|Skvr_pf=f)A3yI^n1IQBqF`2qf@xOqevl~689na`ES+o}r7 z5ow{r9xm7?XgDM(X!ukKD+w?y%wQj~nW@c*l-Axvesg(lE(N%X{Ps_uKCwLCo+y)H z%1{MnZN#x(DB2)t=chPrB3TjS=R9CBvw^M$X&wt3TLlaX;t~?3F+0Fg7EDb|<@4AB z9a;oFo=5}&Oi@fAXF7v{3umZO8-5h(q1w+8^*%Rn^pJ z`wr&aQI*aK3}Io8`;pKtkWF2O8p)OIo`R;ao_oq@ja$1sgJjF%xPg)#Uh0e(+Isu< zS1LGTUdJ8o)HApO`2~7cZm@(=z%rDMVQ>HpZE1rH3{J7Gy#(Z#`1ttl5fP;;tEi0Z zFSe^e9w15qfgCJuYCm<(k_2sM>veBzAU9H`E`hNA(uFoO;lD_`DSYWZQLh{>Y85mw zWYp{(kZViX%fARqT2IWB@oix3$d+KGjb_)U7Z8X857ljaeEiM+bg77^&Ks%$t7>`b zA%GK#-1ZGofV6IfIClrK>&D+2S9Z|ssQm;Ul@`aAu9@`yn#fI^zDLH0|^E;WgmFRko5-WMes}cAmyxE zXlPo8YCt>zF+yp@Z#b*czC@)X>qM2HwipjS>v;yT5s-lFYI_0g95fcloz%8ZRFtA zq=Y03TwoIjmAHQX{>|P)a<&bOZpH+k33#9QEV|L#CiSF$v?ejbL3?9$7V0HzM>%yp zbyeHEdiT;Kyr#x8ZWZh#a$JwUX6AGc%7VId=_xQadzKNK1Pju7{_z5mGqUWw^)uy- zC*?H8fZRDMv>NqT`#8y))Vpf7>fg=_w%B=a4f)-;`MT8n4qy(0@nREvo~yvDECz6m zrknLAflBbDFYy_d#o!|Wm#x!AuMjexK*%zmKwwnLQKkbI-|JP;P*kpJymUCmlP-+3w^aG?Atnxw=DC*^dI z-7j_k{q(8jO}f>xxdNI=Bt4KbHfVbu@Ckbm;+;C6qJoc@P3Kj~E=2mM>3Z-1%rBAf z)sYHR;OAfN(gEVQ>BrO9&!{UYDG3|yxZ{lIENB=yx9-f3vPK6howA8zePow)Dpk12 z4J78c*KK7K?ZCBIMR%g+0gLF^GX1aiF|>x{gaah+li}-lXmbj zJk~D%W!!le6bHJ^XkVcGG4Syl7<%E0nv$6_IS|F5#NMH#w4Qbam2qKV!ELLk9R;MV zas>930SQik)5Uh(`9>ai4E7uTU4<_{ZY?eREaVJYYe4TMyNnb1S*Rc~c$f zUh_EFk%NibJ3g+r{Z|1`u%B?~=X{+3el6YkjczyI1^|eIaZ7E7k=5VJgk34> z_n%Y<&V!7j*DB3hKV@rXXlY(t*6F3iI;R)fyqAaWnWuZKsDiUU8-H&4sx*{>{|6wV za}bS2z`!zsOss)KOy4r`=O7~?C6hRwx8DVPE-MAx_ENMe?NZ}FSGak&4{l~ANI4l} znXLpTg7*UuTZMchpe*JDqv_}*Ng(~-_rOE4!|LjFO* zs?j@H(oAI*L-fcyf##Ww91$LsU^Ylk-)uajEw1X24zPQKglV9KeW#y=*-!!SgY*}l zc>5w@8*CF*0CPcD`1l%k?%dg^*{w%9kT*c&RF`G~B_SKI<`0ljazLTZglw4MbA`}& z+fV`~Y*D0QvOU|N|MBC;9H{e*P(gUSPWN!S0Xf`9A}lCh$fA?GJz0?lHpd+BP>38b z&VeC~5xhEO-WM*;5MT&^{wWwnYgGsyE) z@Z|9 zA+!RJqN}%8uI}o1a-$=h*7EoJ+g)#-WP@3gr@591duRmarfNp zr)()<%1G3xKzUp&b4%RofJ=@?DmR4(jciS>f89t?mhO$^l0<4z0K^eZ*zDO5L9<-A zJ&+pI5(S(ikv3{@i2&@6ItM+JDkcneE#q_Nt-t$8V<0zWz#C^nO^5_D9;kROad@Y( zv1q2y=#6@TwxNSqA7W@e0&JFsG?5{9ntKEGeb!f8jgOF(PaV?f_gu}AOr(O9l#~R8 zSj!ST$%a$ zp~S?+8+cYM-J zh;H_VY^1d5|2qYKWUBwviCsa237<3e?joSFA_dN_>1kc?$mc*ml!Zz1_`Scs8VFtg zdDB7Ai{mou^FB;f{e9{czH18$C-%P!*L?9}akyrgs$%pKH+Nc4#qm-IhsklWqi|Um zDu`N%shqYMt%^f+&!>FSA2EI#)3xE!Z7-M<^}m~^ojx(Ub~*IF^hlnSs&6nM`bEIs zXU@sF6qc$(;qk%3>>uhI*)oZq-4aT_RK%w)3g0IBg~i2jAO`fk5o*;pRc{5`3Q!TaSTKVKpuaxDixuDOVwx?(_`TYi{U7 z-mROOT8q?Tkob7e=HuBWM9xW&KcAcS_ehE)O_GQ zARFzh{5uvPz;HISdr>&yv3h$k6y7#ZYo2}95Ps4yCgI*dd_*M?G4OqyXd!$#2;~V4 z%{xColzMPi;)H=Rj1A{=Mw0>5&+y50PXr# zyGi=Yg%t(B5K`9fiQ@@^vp!5RXZM<5MFi})y*3Zoi7=>pAkK@vB{r2Y5AxT#^rLB3vEwWw(lY4pEs-rlKcJu9dpOD&4Lg zCTH*Nz5bIh(6-(4jHHRDuv*rX6WTGasg*^Y&Y%LK4;O+Wt3owvk=1XtE(PVYhKj&J zpX9cHxmgbl!Q%Y<3$Ve$K3?+27k)rT_Yr#hg{`dw8v!T^ta`OEX35^8)wTelYL5nG zEZ)~XlZ$!s4yw^ii7$KzrWTX!kT0l4oo94gvko!VBKpnoGFFZS607s=Rr5dAbV|9d zYA=+g)(04V-4k8H)%f;vZzEjeH;!YyMf0D$-%~9gkA=GLXzf0|DTANxnln@0msuG` z&S$f?RnmW#nD|q4baX+95#jA0XnVjQ6bt{n*$(8p9egw8ySrJH1UCAJS7MxLkhemApFUWtQc|fpw~+Tjvo4 zHb&`aw#0c3O-*g3&o!iS5^X;zd)a5ze@JXt1Xm9_dC|nmv@YBr-DVMc!JRYx2nic; z`)^e@c*MR+WC!gm|Gba5T?}2~cBtWH+deD{!Fm_*jW=*!0BPVApQNRw1+OhFGjkyL z1|&i$q=0P${*mVC(-aYK?dZC~KM2yB_HF|Y_X%cu zY;mOH8MWR&4=08UWCuz<__w2hu0!4Hd$y|xy^`?d-b-+rKXqCo066;tkSx*#3YCul zY{+1~<_XEpW;8T31aBnI+crR%JpfUVUM=Y7Oo8;^2HV8gY;wnJn=M7yBDw)i%`ZYO zfh23?WPfM88(Mn|L!uZ;zd|g90}~5(mc*jC1gA5xOr3D5ixbCut?@c6X&%~k5qqOd zmRdAlp|Sis>=(h8rG5JJG|Aih;y44^RL#z~UycwGcLcWiMyUny8a=NmPt$@%%btIC zBW&{FW*zzplK*c;366=s0TvAYz=->&nOGJTjnm9vdxV);3eoDm5K9VYhX%JJ2&Anc z-^&D(O}0h9*j?;0(~UzIq)aWQfv*HFc_+v?G1pd4M2KwO(dJwh3&&xF3m0pB?r=n^ zDTHwG?)4d-$YZvbOwJ!~;;E4enB!0B$;qI3V;CgUm&LbKjxK@^82=eA52G||UUc<> zNAs5|&sAxaUmaXCsSJG{B_gq&u)&fjr$`M&b5~qo_w&1dR81? zXV~><1HKJQ^S4H0lcVhE{qY2KlzdeY0| zkFO|3FW3O86jxCDo%fo9Lp2IN>`m#1WH7y6A_W<1ALIJzF}g|Ld9m@K6w0+mF@JZ> zN7T6#$Nw%x@v8GW`>{xDWmyy*)?I$%P2Y|ka}sm+RGp)VC5eF}@&;d#YPNPmLubDX zwo1a=EiTph#AsL4l$Y(bJBR7B$plT7#wFo6IeZ@!(5y^L98r>5N|wynDQx*Who3a9 zAZLI^7CW|WEoDIVI+de00`JZzhO49+`5!V^0!~|A&J0k&T1Ur;M3*=gwZkMyiEg5` z5Cjp+eZxPTKiNiE$7i&X0SSbv@qen(bh(C6DaS3p==k9Ur!$hHwb)rBVcc$_9`V_a z1e!(`SfR&tEnn|1qqvEf3Oz2D@_lZgeP;O3oW_dj5%UMbv8mT2PNaV0tpW#=+nq_S zwbd`*i_o-qhgk7``>S09H9|LxXQ`E?Mzi?$t41U7m>O5M)~UeGIgaqnbTq}!=&v0w zU01yILzFj|KOL!&9Dd_nGkR4=bN3RHhN>&wC4xlR7Q?pQ?Q4XGNLP+y&)GG{lWt}E zCSK%32t#MuY|w>y-m*=HdGn@6<+rsv_~zc89|#rWuqCIc7U(nZzC}%8B1@Xs&8tg!4p9$QddwzNp=E zqQI;i*cdfHyB@*#n)xQ&6Zx8K5UrgJUOibg{WqV~>!8bzJjF1Vsp}MI5dEA_HtQj% zbgyqusWlTtcG~lsp$3*gTMzh`tP`FY;Qx1d1h0bOLh*T|;uyx5#zK77zZ<2Q`kEVZ z7~zxDg(B^~Z~1>f@<@txo7S8+vWJIECqn)Bs&c8terpi9bHs8dOIKsBooi zMcjWN2T$G_nQ0JW>SMxL-@s|?9xET>EXC8`K9@d64WwXz3oGOxN|)voVy~jPLeqff{`b1eimAy6-?A*-K`NVtrkbmM zl))pG4;`c+Xv?lYx$Sq7B+@~r%aI$ez~LaS5_?}_=_zb!6w}?95E&hV%-yk-%XsA? zUbAkU+v)2txsI$uO)dB9`<~%4-`xAXo_H!56=`807WxfE@vjhAy^j8@2VYl_E=Pu! zoDX^p*y5vr`62>2(9mL5!iv&OJl>~c%~9H3he3Rd$j4E@KK1c z8K!wUtCx?)z4KJxwpdKSq}lpt5?e_)61FGzWfVnr{A(banmWcE{+FoF4aOkS7Z?6L zbUSEMSGJogn9+EFF5twCcx$=21s8UJqlgqXO7bx-Egt(uy4c>9{zTmGK_b@T_*;HW zHc*vkEv?U5CtXTIQjQ_W<);0 z!{S8bm2wk;Ux@!%iG~cTfKE+#PK{b>p&!6UY6I~>vlZ_A@hgL?gSsdFGF-npc#(4ieg&k)YLFMVkY)IJ@a##JE5-vB9}F6+{$kON1*%kv0gZH$f1IW! ze!&>`$3f9GaDDhwV(!jn*Ix|l*a}SNuJEr!1+0uJjGU(`JINcy z(2d1@bxh@j9LXwX|HwaEB#XWr7_FF09@?$eqT5%OTZS73~s%NAc&+V>h|lkJ$Fiiwu> z1M)3UIzGP%M<$&bj>pDivme=3rWdnNV)W8*;C-G;T*IG;$UMH@=6Ke)>MVouWaY?b z>Sy==Uh;kEJXA<;x^ulVP)rv-;;x5h{+)wL~plOz>srn)M+%&D(K-L4=_u`rV&yLYkuf~264QTO2 zXh7b;$u6Y!i#!wxYK!m_UEqKKqYYkGR#q;xXcO>qd$wC`>=eHE)FPBV*Vtn-@hKzt z5(0x50t35ddztKGj&bOZ_E7<32*xg&Z{i3B1JN}I-Kdj33z$*gy|b{r25B@9c=bE? z?ukPa!zHVcH8j7yomk+0&pl|s*w|yu$!8*+C;Cz~{f?4}=nc{{$aMM>7Zqm3tHvwQ zHtjJnAvn_ot0(1qXXX9w1@mtY-07*Lpb@85)6t<2NDPG@aP|agsD*JwwAE$X*4t)N~EbFj;kazA$DY*Nzr3izmrI zEbY$i;}g?v&I&VRZ<5qiTdU!YLPJMs6aKJD)L8tXN4eW=qZ`GWEA zGT~%P8IGla+_MjdCikJAT3B9=nzsk}*a+5eI9CljqO%1gd5f7^esnyF6cBdTlMdlb z1c&{?|3%q*$8*`Y|Kq1ETOlh#!w5y$JB6%5gshOgclL@TGO||@N%r2eY(n%N@t^E_YcJdWczjwaW=T%B4vUfxo3S7NOG@q5MP)1%yz zTzB&riJI~Tib(^^jJ3x~L85J(V)b(eP9$C3zMEo-*%g5a0{zF2 z4Sx`dw1b@Oi$}G<_r#z@4ef2a*vT@R#sT#_{zRom!vdE=CH|9_zH!x-1<5cd5l{tZ zT3WYBWYZE1J)l0^^YtlP&B=DVeW*+>_D{(2K(is#Tf8a{QY>$8Z~N6DQs4UbW`l_} z`U3RdXcYe#J=ju>Uhe3(x2uLAkf!bQyw<<$|k-R^tO7BJL^uPtXWomA{1 zmzm#8E@3K-U3|9J#C4>)5kEV~iE9W+I*Qv0+z7|xQ!QOn1Zr^0leMA%*PD-(ZH0M$ zd{vzUnv#p8q+!!{0F>tw6r3!a)V)17D(W*989#3TXJ)Ah?Rd%ECe5`u#_#H*uXqt~ z-MI|BxAJge%?Q^%p0qK$WZ@5omqLZ>?Xo=>Uus@>lqLHI&$PQOiDrnlIMxU~skPW` zZjq}y_pnRh@CDhX?=zp;z(8CdA0LpQ4>4tYjAXoZD;Vm8i(*6%EZ>$Xv2>lxhF>vJqX*Z49{?;uWG%to}4bMN|Wy5BN9uyG8iCzZ_7Hagk`1lxC zL6@JGcN=p!cwAf?08{{lpaVcM@Rte>Y8o?l5nzI!9WK%)tAA_Y9C$K%P$oMaKC?V$mk!yO5Pr zgUAdv*hk2|2a=)QZ$=G8#-^sGXUWY!esDQ&>e9Kc-T_dv2{dOz1%~*5_q6TqtEHr+ zc`4xZe--cl`nc@&kzS*9W$LP3W6&nK&}9{xo8iN=D*la6i?2@^%os69ivP@U6YM(D zo~GwrbG7+-Llogjn9yxEf8%yd=MS1WXw&v-vbJx}K7+zypHMQt-3_QHc9;${jYC;; zZ71y-ZOw7W-dhSiM|(dtozAUWkoTlpB-fv{m?F2nuA{XybxQ>Lq{8eh{k)T#!=zQ( zD7#3!BOZ>iv)^e1p-XQs10A!@FC@r~@ivLak&Np}>rLxZ0rO#>a(CB1nVEb<)R>oU zN)QnX5x>D>xcZI9p+ic#=}q#;6>Qzd)#BS9nf;haP=O@`7$?IoIT!6GBIj7EeKBt9 zzuXHMF6l~Y!GXAR_&$sCzSi2-=Jah3A!lc?_NcTRt@bxJlYeQ)-$wiKuD$2mH~#$B z(T-%;Q~M3Yhl8>!_#-tE$F$WsBocG!)AFzNd1BP3lrMRy5v3uwh3zW{ut|xpp{2CY zHHS8SygxOEN$(!HcP}v{6$AUdX{^pxVYhB;IBn~XqW31>h6)oFIj&&pJ`rBScoGnj z?jdhC7YBxP_)+y&MT>s{9nF|xP?Cz0nkJ=Yi$PJ+P^!wcp z0GIRWu6espKm?BE7Kc5ziLZ$2-5R82%o8TRor}`IFm(qJ@(?D-<#XzY;dfa}duMlZ zo9hDdU|_uLBbFWVf$l{dxuW91G>)~kmyU<`e0$+_8gLlSOc1eQFZA)e=CE$Z$8}f6dtmORPu<$($*NQ^&K*`xz zsN+PL&EpKB=auPL-s&V5$$gYB+K2X=o~jP?4E%^JKldZS(wcwqp=agKH_I>gFq}d} z{dVMyOt2~TufyV~czd~cZB3ztxU=;_&AGN2o_6u0cxizMo}^*v_Uiu4Yu-D~g;^ey zSMAO(FyM*kaO!}MSbGii%{K!de)wXJY9S2(7QNxEPA(^NNY_J5`An~hb%2iy5BCR8K8r3rcaFh#kP``kP zYhFonw&8n42G%4}BEr{LG5tTy`pY7isy(V_l+Ug5P%CNG=QGv8KJ^>KzV_9{4%9>Z z9hHk0i?ZW27rHo~Q{GSP4Z zSW1)Z*+nGiTOp>GjEuUXS?@KAN_!T+0MX6O_g_f9CVuhD4VtFtGSANH^kUtC8D4$j zdcG&AD%2^X?xB9Dy$i|na9`@fiR{Jp6RP{cO|$GQ7Vl{t zZ+6eJ+_bZY3z5Xt3uV2W6|9|W6tBjQ#7P-vamy1Ya-1R$sQ>*1lBwBH@eJP&Yqm=<*{3(ZB#`V;@2Hmj;&_cx@?MD6vkl+jwQqX7|q zqAPzO}3?**X$9J|XI;oy}^Y>Ac zbIy21n`;}FMOObV{+F=ImG_zXk{YP@Vbf3dXM@lWlQ7WcD5VNM$Pha#hU zRPP#35nQ{k(x2>DulCh~g{LNnIUW%+tdk?g9umjFQQlW)6#fp%8-!%4dcSBbh}@cM zjaHRZjp;TU_q$KUEJTjv{Q?5*e>h==jbe$dD@$>ZPEO;9-F7h&o13T=Lf+%U315l{ zdkpjxP2St_dcS|&<#gdQp*MebR*hHJWjk^dy;1C}c!ZkXIESm!*DVm~_LU5=Rxy&W zDAc~ab+R*n!<~Wgu=9sA&W(vjts_UZm0zXM9yeDi>@AAJIK;x~onQG9 z6Vwoq6^{({^)!uS<{**(yTo{{Ohtz?4P}idp3~@=tFsbWoYlrjWzDp)zu3Cgu<>Y5 zTMd5V4h-crm$;$%XZT(e(}oBYF5?-bpDid`R(741SdF@T{o=I*ev->ZiRZaJCY*?4 z`BSlF!!uE^GYwSijq`LO=eQtdky=8QK!8a~fx3C$a*E`IwXN@Ohm$=@El0+@wzl$& z!yFLlS^SYI-^w6zbc4T_=MiA?t3Fzx*IJ7!QLX&<>eXpi!A-7(I-Iu(I9Ao4h0nHS zg=n9n8B=Cp#aFhOw%7jMUmx?~+TR()Z7P)x7v96fiNcLnpFCXAX(+ejk2u@uHghLr zVv>1YeX-f#%2VYcZ+g8CFA%;IrR8ns+sM4zL0PD-aN6#0yz9V03ztUdC5cflE z0%2PL1e217fR>iI{1{h@*bGA~gtE=s7k%!&JTAjJSuz>10N=tlcM(yFr?l|E5~x9@ zDr+r66dmxsnv{g7RDQXK25Soq^>!~VAa*ZETCQ+Z z9}ysETN(iiOb?2Q@2aCK;4<_Uz<8t=Bg$(!77Lg#iDK|fTDMX`_br}DTPbIqskyq1 zhErrFNl}Ycm-5Qj0FRQvuS^y-3Ja;lAu5$N!#~PNT~OKE|AkS3MGeB?-*>@AE#|f! zI%LLW6%{BUHpt=WW7Z+7%LK+Y#QP}eIFQ5Ue+uSd%(t739yphU`)1iLIdqsDO!nDb zhn+K`W}IyIcLdvN18W^vV^}kvEGZYC_8Mu-{V(v;UgsPL<1rw;MDYkJoi=^|*9esA zw!qFN=p!H`oCgH#1(@T|@>v7H3Hv87>41(9t8A63YiDzqOLg!I-ZFeuwGETQ+)Ijghkj>JFzU(*lTzAT~D` zUtC?41LrUxx3z9o{rvg!?jizY(Zia>K%6{GuMEY-q3htT%4}B?x$CljX-Yr*pEolA zbDACnjZU3dwh=N7E}I(bA_wdO6dweMWPQrod8YLFkqro6kCUP&_@n;x0;%&q$@f*L zazOSi)#WOl%Be$#Q%6fG{`<ssdsKZDELk_y6NZ0iUuG-p2GKeMKTSJVn zgS|}L;-gI0J;U9EZjs~Qe-KodSU^yPB`opNxo%t{^zz(a)f*it-mjvgq>Pz2J0aP{ zILjZga2+egB(#*|Hl}y0b|tn9+ZHU-bc`%_w*SXQ*PS=Q65YB>8RuLgc4fqA%t>#u zhcELEtUEr|>l_1_7?H$J`T6!M1DzG3Flu#mVn8+ttk-R^KY}uN zofUd4#j+@;8#rqNp%t*) zIbgn&tZV+m<7(fe1s}kK#0@RE(SLGo^QXnENJR`KfsD*lQi+bXYWs5VGQms3bUB{^ z0(vd!?iZevb#}f(!RB!& z?6VFWFS5JgzjOhw9%$-#Kt#|7@iY+S?!0IUeD2|dCEfGs@rk2d!<$RPNma025W|bJ zXnbnSu~-gPwPlX%V;Se^BoaeeU@6>1-V^83p7y+OJ=*)~9{Ydg#nd%U6D|P&fNp~Z z1N}Cz(jNv>)&qUa_-6_lki+z0tvOFIylOv%_F_teZLf_+gq*0+bS}hybN5PQi>lkYXnj2eFD<7aRx^hkqQBx=7R*qgz5c{9B~`?XtbV4Y`;)Fi+(Yy7Ax8=kWXcdQYnb%*XMh zLSK502#ASOIgonUeJ@TA=B^%nd{TGeSrP8T5I-D(Gg62+W^|h9mBX@-r6*;R_SYta zAKds^YIpbv`|%tT|1*pwogw;tOf|fZSROxSG?eq!Gc)S@pW+KRt+nl~;4w&gP$ZLe z%ZGZ()T1G}L!Z3lQmi&mnQr?Ej7tY4(xhAKG)zqV^44&Zginl~bZc2h6eK)T~ z<+Mgnv6w<`6j?q=N^khj#aXW$ENDp4WMqcpieSzjzkQxLW}R*~#jv809&GBBIq7pCeS5FaF*>OO9&u&o6Al1fke$8L2ZGnXr`{nD z|2cp6r`RKqkoJoOuO(*c9nD-+nDEpsoJV)P!uJxAR9b{Q#`YfKg(z#6K0ols16 z>yMk?XJfwI(@LALS=111YKsm~%@$)hSk5r7E`dQ)i<0@1D`<+AGRVG$q4GahwB2P}tw3Eps3*kKsRrRU=}DGm)FG8`>Nt;h;tT*+}F#Nj1Kec_GyKT2NM-0Gn- zG9pRGXdb=i#Zp0FInMY1qSUpm-p|8M7GQIqam!ci*V=i!Ln(&DzqjOP!in=f#XmkD z#B^&;bnOlpg4k-*u)XfqTWv`D;tfdx5BKIWqxK@mr`gOSnZho^wAel*!t{b%nJ;Z( z6lXIj%sf<-$RAJK{;v6Gl*Qm>iNKRDw)t&KKgX3ZP%_+y5gjqQxu}O^!a=VS#PJ>* z#&qkJ=DB=%X;jViman)U3A%rQJyy{G5_iXLTVqdQh;y45CZY?ziSWXj+~#w&VSM3u zV!yAT=~gfe-%9oNJjA+qsc`-X)*bzW=kJgMPNeEoEVoZN)A;(R!K0xBh_#JWGSPE# z-oHxyz2`a*Br~)Nc|D3QA)Xqd(BInKE`E~Z zmkoJ~EjYkANEnZwoeqbmLA+tfZyKtlr2Or;?@DX>V5734NoX`@chm|#Y(4`R&dkHY zyVy=&e%7n1Uv%+;`!C)#M!YOHCd_3xL|BsIOEJXMc0AE^U$;}Kqt}U^dC~9|TqlMX zn0#i65SAW4zA0$R&QfqO&^-cqJeTmA(0dgh5l#>)b8$UW;)ijkcXgPd|8*>>i$Cf~RkE-xEth6yxA;6o+FMQTd1+i-eyLEfUTeku z`Z{v5qnLdInjWkM8(t@ATob$fJ;mW;*axvy0QeXe`2-}YUm;{6fi1SSuxf1^cPZYN zU~krRQttj1_s6r|xpA_q7PWH~_S-vHPN|2Egej>M$jhNujH^DVX&Sf0=191rT)#Ax z9?ugq+eKY>oSYibQzwC239Ea?zLHI$K>Qd-7^-6d8p{r;USmben;cNZJI!> zsVMExuR8R|c#f9(Bi8Ln>6};T?_V94u%G2pAR;Q<4<~Cz_L4mf?!xruG~N8JygJFq zm$Hjir8I=ua+Smt!0i^BuQvL~A3;Y0c80`(->&Lv1jt5~Ki%1w7h<8b6wh0H6*X0U za9&@x^3>M;$jnl@sQlsI_n%dXxPdGcP$EZ{2Nnt{rk+uH^{N$zSrGRJf9w-a78#j+ znpF@<42OjxqdpT0Ie0@o2w)sIh3Bt9J%7Sy?8ipL|iw}`3&T!Pd!GRO{XD*SDsLoES%z7=Fi&uy|zpRS1; z8>jc>{X&wq>@ic^s(7ONjYdTpjk!el&~R{;`ej{DMy+DI0^o$c3FV*1y=m^2+V;mA zW;*|_pqw(*?6|$F_-MMQTHVK4xhED3JGQ*VEP5SnKwlu{?7~FXXt52=-{vQfPo(PAbhXGM3qr$I?eW7u_2z8g~m3gx)^EX1MoQ_S-!vpB# zfuRA70o;iS0eDyVlr`{A+&ZG~jA=ig`xwuakdWf)w3Z)?O(8*MU#fm|6?)EbrIYY^PwJ+i4!U#G;kR3T^(T8Y%paDk za2X`E7a{Hlb-ww??t&Mts0AB_m>+D)LBNd+j`!p;NDr;9Q7Z9e4L&4~=8Zz9vk?Ko zVBIhkBU53OP$ImLA#5TQ0ZnW#S+bi7s%WI8T4}hOI zZ`-x@*0e{$HP29X;Zl!Dy!2%YuD+nlMf}-34=vpe?v9NY*LbCLUVz3bSSml4ud>n< zM}^~N^bT}CR>zVg3q-8J2ndCGctZ4ywltVK>E12qapD}egR*v|87|iRDs#z#HR=m9 zpDT~D$=w|@?7BOO77BG=)rvHop=wRcOCJ0nlx{>7Xh@Q1;Ja~XK5C>=^F<#0%JkXf zwV&YTf+vYQ1Jz$ub#*ov=^x%CXH!iBC}v+c*Pu1rq|46Qy5%$NCdx&&d^j_Cl-)B+sn0kFrZL_LIneO|q>1?K(NZrsp-)I~!l z2&5-K`z%a4L~m2#0c9-g`)0jC4kcczA1A+*V;|)@^-ol(y~6jqUmPF(?Ksh1f%x2Y zbMMEUtH`)FM%p)Z-&fQAa`V@7D+)UvtcL+osI)+3H(-|Kp%OA!lcnHnrfmGqQXoz# z=QEeq%l}F8yzvnbYs;IP-vUSI7q6Q^z8nWOWl<3IL(16#GD~U5moE?zr9-)ZDrx`! z3A3)Ycz9D#QU;c0vAy%ASsBQ;SsRi2T(~xAL7PQMC_ZH+60pOL@rI-I+elGR`vbZf zo?MOouPVaUUb6TUd6mvvs3$J84TK`5Ty31X!8E-kEykH>kU{03=BUr`APM{q{8480f@?6<2TaM8p64EqCUash@_#dwwTN%^ zz6yr{LVF9^RvaxBg>(4{#gr0QQ;n`*V*|ET?^g2v0 z0KP#ta|g=8v{(96&7l|+b$NYhP858FQN1g`D-+E0g5m##S=t82?>WbkwTREjhBp>5 zz5$&}z}_#a$8s`#u;}x2dVri~QKs(IHT&CGS-jAsDFo<8HRJp2gdAoVqI*%jHJ2Th zjDHXQic-JiBI$hr@%8Zu3<{b75n|0^F=UJ=g=u$0csM1$?KPlNq=WLve6-{p=zrge z#2SZ?#p-6F5#alll}U+=L)pB37rQXE?I}Wi?cr1mI^F~0KTI6lgC8-<%8Yj|wtc-S z5s2Qvg(cUBxmIbVh_zO@H~w;C2G0BTZ*PP?!K;JfMS=y@>tNl*0%ji*2S+m*G!Ai~ zjkeuckOxx756~)r1E40y!a&yt7K{&UzQYO2FQ**IM{GXzvgAC$RI?-{d2++f#l{Y! z-7Rcoz>tFwzTWz%GLp+j%rjPMZUM(ZK0>kc-=`J{++k{DG4RRaj`}>rBhqlu(rJ@( zmH4%RW0?w8O6!ws&vx&(ICR4KIaG6qm#>Nzk-EZaJ-$P+8Mu=E z79FBZJsNcagdIa%Le{u(KYn8UXn(k*o35HZyF6ke=lpIVlpv zYeh(GIDx`uQDl;8K}_q#rbc&Y55ek{vyoR(@-r`x+iW|j>=$q}?22kduU)TCmk&+2 zjQWjR7-G~q17{z{m$L}1Bhda<^6^HQ)i0X%_tH)Fi0r2)CGqGWj=iL zygyygnzu|;RANpaLKNdyDR=-r`s83vf9OX4&&4Gc#TGX(#E5SaBaURf zH#RM8+b;bSIpO+n`Q}YNJ8Xr778MUzhXX)9PPRA>CXy-K66)ryS?l#Nrew=2xJm;AVtukB0h;O9K4q2XyyC^E7uaTLKk>ZsqR2(b7YgOk zza9Saw@@nyD%5&{{`mK+UyaPl5!B@6;3}<V~HQ*hAF!6 zxjbeTS8Qy;MSRe^FmA(;EuY`)1rwejRf;1t4QwCUY-hUE95TzXwD}}zu^ZSvb)Kh) zukZAacT%E|Et%hfT&Bd84fxlo6>SMw&ICV?393u0oU6XpyR#`f@dAA>tuDtdGXXg6w1&sNUWn&|&%T#l7HRlbFH7xg2><`=rU zNUmh>V!261y_|CRpo}cps>iSEzH0d+wSn`c&%_J^0~1@sRvU&j>Kld_hZb`6i*5Ms zJt~B8{RT`ZJ1;LgF6Oh^voOVd_%d5p*P9#L2dn(~B`xf)W6wes^mFN_RCnlggrKcS z`US-KA7N!2Bym3K|HZ7-b7|~0?R&JVjPAfp8(kSwAAQM)U>d7^ypsPwZq>T?=)8DB zhkwV#WKw%{eUcbl2t=;H3N=KnklGA9q>|v?ZSp~pMgt^JIuuxpdjbxP;P$Fj}T7r$nIE4Stw!pRJ~Q#ftAiE zjCu~JNZ?mlj(NR*|IfK|D7;XZyM@MRN>^8(G7As#@=R14U}lLqZh*6_PxQ7W&9#yOHanBIMV?#%tK<+A~jSn4*8TwKdG?dKk4Hp=; zv0yO^tbTbwIuqFQ!NB#j(d+JcP%uRQ+plk7d!_Vamat6TSfO`h*Yv4JjJg{#mDOcs zaT#puQLGXNG0olC!I}O$f~MzH|K6JsK;+8xQKRGf%j@euzs2r7i}MDuUP6z-T>BK2 zkc~XeZsn((3kdXEf;$nQ?gUL5aB+XqE&8XdbiLGt0d9moWlL450{&%P{8hsfO9&Vs z5r!}5gfT#&U_vjJXd&6tJf8R<_o40VLuvV93s@dABB-rYRX-9&CW!6-3~AsiUgS&N(Kl+WMhG%R zCF4B5yQ5x-=w`3N`;9;5d4&4Jd-Aczd>4G4O*jOHrA3fV4PD>l&;fzA8cV7NIg^a%8 z>eT0@Y$KQz{S4_vB|hP|PH)hm;L#0}U${Z!zGZy47xFE%Z=E9g`l$ObD#9ytTsV0R zdV%GF8uK&^gJs?=V>mlqLVX)OeZ|d%sSmn5yyAKOPMI$kDDYTHB}{ z<)QED8mS`#MgoIhr{wc@DI=Lm_93E@tgmV;J&E3zpGz*F)&_?xQf5K(y!^#z)m^JI z9n?vP?W2IRH6VJ1oU9K7bN4r(=rIN%s3B!lBSJ?na&Wc~rXK1gwmx{OATGrvD^r_a5qkhu*lzO}Q8^g$&u9*|jzINjskZmvK_GEgPK-95$8%a-(#j zIlEYZn?Ap5(tZv5^yy_@%Yw;#R?YZzG$_Wkk9jgGLpp`~7(D4wDkM>Gh-%v$`p0wV zA((=o+Qfk3TLfhg14PG;@3e0NdJ8j3P8aFRz!|l9A&v1t%0t4n{9ovJ47NQW@l-F` zb)5(l!meBOk{GbsV7Zbm4F_B#`Xxj_KwxHV&Ek~h$tTIz&0u~Ijl%RGh(s~tpdH#% zl<_jiX`|lqQxJdF7T($ggLGd)SLwW)276#nGk17Ll z1IYLd_`I%FQ;UIo9-i7KI$CY7HVmEdKPX*{-r8^X>({5Br4v+?@hS0GNT{b9v|2ox zURn%T@4XII-s1O6Uc|cz`Q3^xZy`fH>-+XGV=kq|zBtJ)>rc%2LOP}QWL+(=uWpu$ zX+)b@M7aIaxD(Yc0j#eoPtQ0=?s+-cIXPv3W*m&((~c|t+t>n3G14{lji;d}LSf;K z@G;@PDOxPYOug6-K*z8q~rixC5-0Xps3ADhkeeCy-*?3yDS^u#jj@}uujGo9a& z;isVYll~p^quWw{$PY8)2hN{CRXf|@UVQB&z?_HoCWkZY2++%wU%X%hOQp^{mQt|v z*jLd4CUX)PYPZOTsqL<1S*YdwTEHs$@Z(+aT32F8yDyOV;^D}8%gZ^NJxWg1hw(_4 zCk>rzTg@O8RGna;5Du2hm%!#7HV(=e1R8pPxe@wpuwFqqA+KyfvzvL~$D(Q!&~qdb z_^}1L_2>^nB70O^R5KK&rqLcUtxi}}y1@bE*P;h|$?9YvrXZp53+ML&juE2}ps){r z7Q^&LHDZ?8eZ-kYIF4;uo&LmbfWV_dy)<2A??SGLa1&VX?p1#OqW#;p-tgBrji zrWL5=CRwFd>6*S{RX7dOjoQ}ycufdlSBrF+eXApse94|c?89uId2J(k0>2Jt_oQt8 z+s}{sZdqGAU5U(4*P*XVhJsE_E~j=C2o%Aj;U!S<`2o~mN^-mCaYp5FN}ev)gPuR8 z`h^54;!>fe#a_a%<~-=28l?f%s2`^n#hwhm&kcL8+mBxMuoVX9BnG%N}+I(TG{WLT(F+fppS1u7!^-o`qFJ4goI- z|6ByV%!b=j!n7=j4UlTA%Fnydq4mO2=VG$pVs3VyJJ= zV)~V|m29nRdQ)p6=kjVdnBZAZ#YoQaO5cb}z(Woo9=$?#5#F$fZQ3Mf z&omcl$I`lRV|P$Q6Lb$PEU5QbT+S$7>lPi&p(@iBwqOI{72{hFl6F_<4D9Hom2_>m z?zY>F%W7ury$#Wk?#pc&yQ8Ga9qJ~NhR);unIgP*(i z_WttcWwkjWh|tddba0P29@>ruH)TJ$PZNCcG7eBH&UWFdi8S_z&6=g}!KYHl!55`V zjI7jXAiCj4aS=rwmOo{;*{su=PbM7?@=~TaTKof<(4js&v86_I=b9Sc4KYKr&l0H= zEwFT-x7&a|%oc?;dD=R@h5~rA=&3=xkhCbg)Qx9x9>-Yy^B4L&!foc}*F(~%m8ifh z81!}p4ix@ZAmy6VfdmKmt&x#qwKX6O43rBjqZaRB^FI){hUcZ|K5G4?OMH3%rY$O} zI5!<6dK3vMGo149=a6Gps$u&;=i-sC`dKl%`h<_ZX3;u+X_{LsMK(R!_dkQI7eymH zY40ID+U=v`)QUSE%zcVaOFQ-ZX*TC45(%L?crE_jWs#)NQ3cjl+Ux7m1+o{BE3Ec; zHtT-EKs}zb>K4s~qr}h8Z>+lHEk0PH2$%F1*aJ}t%I#bxhbjveeJ_b(mPh@@9YO&I zUCvY~3TH38iRz!t=xNrA$@;X{rxqHzogdvZGabDjV0%6S7tcq?M?A~RW36c-nVxD6 z=w2(R;rZN122@d$MP^v-vi)scue}5CSv)_2C9H|sCi{25eRfPLF~zEZ2W)3gE)k8{ z_haX_7I8LEpjL45XDU8s=KbRth^PFV5m@YByFI7g;#%Z22Wk)wJzXL2k~?NJN3pri zpMF8uHBFgpqS@QW%{0;LQl%`o$={JMcE?w3+eMv_2QX#cV*eG)Uwt1f5jCqfbh@0=>ziE}J$*pO$LqfD-aAjEfi2l~_MY!=dvUbm7%2 zU`b{XO5a^=2I8)3FnlrOOoyjZypL>RZ#TD92_1Otj;f}+{RIp^O~Dy)S?3oW=d<^J zZ+KHt!qjwK3nY&jAP{~yyzC)#xX$uBEn$lB`!fHg0Pol(H1>lC7M;IH5c4x_5B1R; zO`M*da=T93`64vtH*c)>%-5eQqkaI)dY^cw#vJG?AI0pw54@yk&vw6P5=(61_>v2x z)Sq-SpHpd+Vi!_TaCRk{3Gf-W_n;;4g$V-vY>o!aRtwP}7sZ&| zJ5a@T4-}F@7kjcPba;)771OAz?#H`wJ9$`;d9t&-Avf)pHB@s2-(G6{wrT!jNR0;# zW!}*h%-6`atDSGd;a?Vwx}%MQ|6)=fo#(#R`t!m^E*U!%K+z=}+$@ zi7rNFXgnf&68Uk(i)<#y?A*LfvYSy9&Lpzuxz*_wiyf_p^;f?Kc~@5M@v*WR_X{973AcCcQXn; zbJsFcvZ}pSNEBm<5#%+>eVUU{@_`zh#Dh^Lzb-ol=^a+>0S_%~7fJi$x9^804(J>F z3}^{O^`T$Q1BmBtH}vF{1ng%RO!mMitYV`a3jJx;(+&JFf_iqjwL~j}jXqenZ5Tlt z*x`_N;aMQMA=nAKQ47U#E*8&p%4&2@YTk?6dRj{&ra7~=LzLW}LLA0P9!OnR3Ehm) z9)sTv&R;+YtTN;Qof$4NQC2XLA8MEly+2%Fy~}z_C&_LwAJdC}{MRkN78}$CIuF+C zN{Fg*F~J9~2fPE$0~+uMM|B85d9EcWGjL^y{X!;WSiVEZejiLTzd&Qp*UcaPztEtr zyn)+y?u0?lvK4?iFj2`7Kvu)Y5AjH3hldJ+|M=p}f)K}aK~5>j%+B?4+?c|xFETpE;ySr7XTMoPsBrPrM)k{$Et~kz z!qbYw2jXSjC|2TI1(wg|iEmc_9IfQ;fK3^yH<6n=eYPT`8_?{^9vX4p7$1a$x9*?h z2U_Am^27B(Crpsw#-&7Uck<@6YbDS;N6mDw>+Ee>dT*hQsQpXs6$*C7eCe>&+SeqE zQrQ({LD!9<^AF;G2+?P}+q?2Jk5p8Y4o_?{+}&kj(S>GF$)Fz?XU==WMgZ?Lg0JMj zaQ+)~a-;@Y3xO$e95`G}ZET#kEJBabSeSn6KTAFKI*R|%{|)_bmltT-6hD9ed_a@0 zBk^synVOR4xv*5u^wP>($aKP;)_epUZ0D$nnQ|k`t5Nymu=)Z$wr9(#XN}%+!pNVR z{6sL3FOhNgSo=hHK&fIxGHKuL&+>V@%*PtViwwHVX*r0d>urp#`e+M!9~FR?4(PqO z$xqNw?F>PHXLL{jEu?@LyuZpzNpQFqdX;y$Dq+I?@0^&v zrBvF>i1?Zlwn(afyt$TI?6IaShD}-CRAl&%)mj~VAcb>DjV7)JPe^9RqKk4smEHDmw1RYO8V4hhWU4LXe+8$lV ztHmH$GVE^2n|e;wdHhx&&)f;O>)KY9l=29n==J*KJQ^rWANOU1mPn2jo82~2*h|;? z^1h^Vy-o=?jGbatQ|#@q#o6RUb@S7w3$C>1tDS~ywp%{l%*^xb8)({Je%f9Pjk*_- z(bH$v5*vQW#qD$Whjeb$W=)xSzKH^x^lx_GOjv1*Zhf}DO)k3??(5hjiHlWLN_vT8 z;84Gv4ma+a4d?L6@!oectEB3I5m7pdgL#($f3?f!e&sIo75eLhw^tC zv^beh-V7=uWGh=+-m%}N4gNny9Us46d4l$MuA*>fNFhl}AfBCbvc=Cze0)bw|J zrjluBcnMgmoWtF7#5Z@orx6eBAE{e%o}=FszxJJQMl^jjZ@<7Cv%!h{P@n5*;dqa6o0E(WUs^r7;Sx=fZo|kkBKFFrXGS-oH}wYIT_A8_ zj8E-1I{e*5$GE$FuM!Dc@#$QspquhS1Et>PQtEBv$Bm`R(*l(w^vOT66(y0uFOw0= zinfhg))U9c;=vr_Xmys8>q54n^*J8N718o<GnnDwtVOz%afM37_CV6iHMlyt@S?M2p}S)6o(~lqGB@I=e(C#-r-y+ zdW^SY{E+H1=H`m!om(}YG#_fPdCcWPlNDYJN(z*&e|sVt#s-7Ctwk)1>3M)*@UkII zX7lVgM;|lS4?U&jDVsG#p0KB&XX_c#eBx%icSx!rfY28&wWo@T*7t;HqDj|WkM%ah zaM|p_0V`teLi4a250vW|Ck;y*gUF{Aap>I;HZ|?TiU8tDIwq!Ec?(AEZ*}T*muQl& z{+?%!65b^UJ%4dIYtQ0(Ld4uqYT16*L~u+=XzCSwv_TsWS*csse~JE@YkLq!hb0yU zPOq~S+_-g-EbkND4OiC2S6;1-?rW%ApDFGw$spaeCEKVzG&(zZ;dxebcD3}mKoh@2 z%!tKp;S75H@%N{?eKsD9@5E7$ck>>($QdP%xmUiv-lBr>#L01V?pGhq8f)Q=!EUDO zL3bVCOjt1J(ChrMufEx(CCAZKSO0`HzmU-$w{YVeC!?{XEPC*>mLzom4AXx8==AQO zR;tVJaE@%#H}##~S!LuY^8+xA5Z&Z_keeZ&dSo!6ht#*7%Q@g`0%&~bS z{Ah5!O368lBHDLSWRj!fjajFwll#LDRvamr+N6q4fKI2>+Dbw;hvzz4QZLV2Zb|R7 zKrlB~7DDAVOE`GPi;?2#&Owb!au?G}9W_;}5Tn)d?yW`#u4azf^^^U&s{@4NZ=TF* zkO|b^Z?Yo;Y37SDK1x0kk3-p^0xXf^j{b>nRJXJV8ehk1Bn8HP4c#bwN;y`_62h)t zT~=*gd%VjYBGM7MKbu;wUd$np{zR0t8$M(vUzUtpTD)>=Uw6i9O0%tyOJ?HzJb60%5hb$)?J{{k~)2~^i7NyeEgb& zg8gvl?^fcGhJH<&A^bx;r#S4%XtDBkY|n+D_T*N~{9@aM3OBb&nUOrk0^W8h+O*p* z;$7Vn``K*WVb_y}u_+0E*ECciRky|<_UKJasFD^?Gc1oj4A^*g>vvDmIaN)>TlTYc z%R%*O#=*;2Vx**p+V?V?Cr9_xepo0>MdQdGw0;*J;ztv^;C;eUbUC*md}2)WV;_OqY1Wmt5sw;qlwrwZEAEcRPo+@i|;KswkA-c zICqqRFf+c;>t}KsntD6K8B^;ZqUXhG2Qi$No381oi(oCX)b)w z#wpuk`9LHr{XW`b$tH;!!G`@83vVf3ix1xC-R|pV8a~?XY>GC*IH&Cq)UAC}yvFoQ zaLUYaS13qNbit!DByldrz{7S{(cUdt`Vwp(#+i9j?F5&zfkPFlr$y|R?LW1&p4nXU z{%Ewb^xOgSU@_k{)}vaT#s$M{HpFge!B_XV?KCdk`9YirYIryw+Zk-Q9lD^OVzdg! zx2LDZJ3H)WVJ4b)PF3@6D=DDgSGV=t{FQe?)Q@A20CqjMam!}q5$h|Jv|0+Og4Z{9)Q=lZERCfQGwffcw!=oouw#p zc#VaHh3?xk$+Vm}G|qi=h>0) zJu|C3U+jA$pDmb|X~!S0Pi!ObhXe|q$7VrlO~)?sh}-d6DN0KOFyU;;sO9VjHElgG zA4ll|zkK;3(7w2^z-2op4KzrXolf525(^sWbCCvIN$B1A!!K7%GOhwPx@NcK=RaaB zS`wAN;)1fNW@>V?YkGtbDHB91v?QZ>kbU{dYUAkQQSa=DNreh+b7*1Ra7FRZNdCW3%&Iw--QzfrTauGp4l!xlhM=e=~JINp@9*uLp0^*%h02o7>t$pB_*@E&9UB^s`L*{gX#*UG@9H6a9;~7-%;`XT;tUNz9`!~?pqX$TF0xKQKqR9YSmM9zpns@;hw=5{P zj=Hz6ual$R_Sa%keyxG%nj((1mB|7z=kvzvZhSJ#-}iXP4zZM1jTLO{bb>U8n2Elo zlo?We|7gS<#5`rNVpX@-C~nn3sJ}$>>vF(&{TtiC?2(GnB?V)3=9~K!T}P)g+oD#V z+{bs~Y?F1$tm{`V$}N(k<;8Wr%-Sv2qjVVYFg`pQ$wKu`h6d{UiwYSLfONx{10nV- z-_W%>C8<cl$B!5TcM46!4h8fKRZ8v09-TR*r%uT8<Ntxc`xmI$Pc)V!-@P){0qN}6qI~iCk2!KQ5HKb~-x}1r2)`<- z6c!dvJYKKNop1nOXZ2D`t?UvgcYv>L#+Cb$)T`4m0aPP4%M}xZ4gNfkxyyKX)ErlMaKG%B0(=mvB-CDjSV!Gb zhe9(eRF4c^zC;x!K(MR75Y50S|14mNXl&xQ>I4!+n~Y$FuyMFF)J{Gp^J*M zaP9>!1P&JiXiL6Zvo~Uv+iyYxtr;qv2CH|&;4cDs3*)uCQ2zs!rlJ1%of-AP3Vas) z|Lb4X#xFGW5i}xo%l2{4TX(ZRee!``NK^I#X5;ivbFS=`^Z(jN$`IKXvRc~teO{V% zPG8E1RhRWxNeNxdz#y&l1Kd*j{)fm&l8YB#I8vqOQK^R1s%3ut8Z_G+dWnq8uk$5$ zWx`5S-|{uEKFX~F#V}n!3U`j%?rv(n3H|CW~9ASy8|W zQk)=F!~!~syj=q0U0$a`MQCM+Y&sY$R;S}lDpUh(t6Hm0UJF(vfKWh~C}J&!Kw(9p zl1b}Efc>u@_I2OaFHM`?+~%Hpo^$Rw&-42^IDA5<8_U&|z8z`rzuxq8mcJwGk2pS( zX<@kkX(?$aRPwO@%1MLgLi4oSf>-0rgIMroV?ffkGRNEmZ?|YByPBu$DXSq~>Kf-& zME{WV>%WIc9=4t48k4`iMHBwVh3a*46U;@Ww=EFH#>S~2nYAe=r>9-N8fHY^Nkt_H z9vYuS72?O(ula;j%`2SB?tg=A%((tPbtJ35`73TT)B4}PPgqf<7UEq*p)jD#Y# z8eZhH!Cg;nlZ6vRra>(TiHz(eZ7DCMjmPmZ2M+O_Q4&Kgyju)M$EoWGz)3(o!N2n@*uX!E#Topx7%8kB^?EH-6{6@1;SNqtBt6$L;&$!WNGvwsNofV{ z{`=UTH=|PDNxDi(jyDDW6|}~R`gUUdN-8$#hSPu)Ym0AG394>Z*mn0!=yB04H+@L2wA7thh8sJT1d=+N~$OM5)BXKYGKOSAP8pnc7Dj0}R*9s#Qz^Z|8EO=yiN9$Jo^ z%^*VdIU@%}r&^db^pk%vEE;31X z&Kf&!tn-l+1za0cvxY%H@@!IeKlg+n|FPGJBSPBNQ$P(4$^{v7Eiavy^^J;EudygC zin>^zs>o>BX3VfalK%z_W`-^0+F{8fY!X2Qk1G1>gbK0@2}b5W%j<^uZZE#yecE&t z52AHpSybKN-~qS?X3prVwjcFgyt%p_bbI*YL2x_BcPKlM*OeU4mo&`seTH(V=e_qrv{xOTy5H!~DWpuwG9<*r~nvN>ZP3Uxi zsQA!owWN?o4QawaoDLpyfEDFfk{alVe>J$3{F10KS0I%VL{#xRWJ|*UHbf3AyQh- zqY{J1yV=~CTa1|t|JBp>D!~{|^?`?AE)jVl&ta$sJ?~<4b+}(Zz(xm$Vr1!%0}zaW z@JrOz)>cGUj=*{HVGTz}M?E+KI5TCy{u1{y1VeZEHrY`~yeR9vUM+jNyO#qA|2u&8 z$p5dvoD=Z;8lVFJTPuK$$;#65-g{KfxP-@uR#=YMd^BSU-eEE!^z1F|KN6(Hn;}8 zvQVK=SXo;iOiWCyK<8iQ!M~1ajt41MLPTet)c+>7LkuEA8uIa?bUGbe66*@bTWtD} zk|RWdSI0Py2G1f37LSbdq0G+1h4HlHf`__ +Kalman Filter with Speed Scale Factor Correction +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + + +.. image:: ekf_with_velocity_correction_1_0.png + :width: 600px + +This is a Extended kalman filter (EKF) localization with velocity correction. + +This is for correcting the vehicle speed measured with scale factor errors due to factors such as wheel wear. + +Code: `PythonRobotics/extended_kalman_ekf_with_velocity_correctionfilter.py +AtsushiSakai/PythonRobotics `__ + Filter design ~~~~~~~~~~~~~ +Position Estimation Kalman Filter +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + In this simulation, the robot has a state vector includes 4 states at time :math:`t`. @@ -61,9 +82,28 @@ In the code, “observation” function generates the input and observation vector with noise `code `__ +Kalman Filter with Speed Scale Factor Correction +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +In this simulation, the robot has a state vector includes 5 states at +time :math:`t`. + +.. math:: \textbf{x}_t=[x_t, y_t, \phi_t, v_t, s_t] + +x, y are a 2D x-y position, :math:`\phi` is orientation, v is +velocity, and s is a scale factor of velocity. + +In the code, “xEst” means the state vector. +`code `__ + +The rest is the same as the Position Estimation Kalman Filter. + Motion Model ~~~~~~~~~~~~ +Position Estimation Kalman Filter +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + The robot model is .. math:: \dot{x} = v \cos(\phi) @@ -97,9 +137,50 @@ Its Jacobian matrix is :math:`\begin{equation*} 　= \begin{bmatrix} 1& 0 & -v \sin(\phi) \Delta t & \cos(\phi) \Delta t\\ 0 & 1 & v \cos(\phi) \Delta t & \sin(\phi) \Delta t\\ 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \end{bmatrix} \end{equation*}` + +Kalman Filter with Speed Scale Factor Correction +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The robot model is + +.. math:: \dot{x} = v s \cos(\phi) + +.. math:: \dot{y} = v s \sin(\phi) + +.. math:: \dot{\phi} = \omega + +So, the motion model is + +.. math:: \textbf{x}_{t+1} = f(\textbf{x}_t, \textbf{u}_t) = F\textbf{x}_t+B\textbf{u}_t + +where + +:math:`\begin{equation*} F= \begin{bmatrix} 1 & 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ 0 & 0 & 1 & 0 & 0\\ 0 & 0 & 0 & 0 & 0\\ 0 & 0 & 0 & 0 & 1\\ \end{bmatrix} \end{equation*}` + +:math:`\begin{equation*} B= \begin{bmatrix} cos(\phi) \Delta t s & 0\\ sin(\phi) \Delta t s & 0\\ 0 & \Delta t\\ 1 & 0\\ 0 & 0\\ \end{bmatrix} \end{equation*}` + +:math:`\Delta t` is a time interval. + +This is implemented at +`code `__ + +The motion function is that + +:math:`\begin{equation*} \begin{bmatrix} x' \\ y' \\ w' \\ v' \end{bmatrix} = f(\textbf{x}, \textbf{u}) = \begin{bmatrix} x + v s \cos(\phi)\Delta t \\ y + v s \sin(\phi)\Delta t \\ \phi + \omega \Delta t \\ v \end{bmatrix} \end{equation*}` + +Its Jacobian matrix is + +:math:`\begin{equation*} J_f = \begin{bmatrix} \frac{\partial x'}{\partial x}& \frac{\partial x'}{\partial y} & \frac{\partial x'}{\partial \phi} & \frac{\partial x'}{\partial v} & \frac{\partial x'}{\partial s}\\ \frac{\partial y'}{\partial x}& \frac{\partial y'}{\partial y} & \frac{\partial y'}{\partial \phi} & \frac{\partial y'}{\partial v} & \frac{\partial y'}{\partial s}\\ \frac{\partial \phi'}{\partial x}& \frac{\partial \phi'}{\partial y} & \frac{\partial \phi'}{\partial \phi} & \frac{\partial \phi'}{\partial v} & \frac{\partial \phi'}{\partial s}\\ \frac{\partial v'}{\partial x}& \frac{\partial v'}{\partial y} & \frac{\partial v'}{\partial \phi} & \frac{\partial v'}{\partial v} & \frac{\partial v'}{\partial s} \\ \frac{\partial s'}{\partial x}& \frac{\partial s'}{\partial y} & \frac{\partial s'}{\partial \phi} & \frac{\partial s'}{\partial v} & \frac{\partial s'}{\partial s} \end{bmatrix} \end{equation*}` + +:math:`\begin{equation*} 　= \begin{bmatrix} 1& 0 & -v s \sin(\phi) \Delta t & s \cos(\phi) \Delta t & \cos(\phi) v \Delta t\\ 0 & 1 & v s \cos(\phi) \Delta t & s \sin(\phi) \Delta t & v \sin(\phi) \Delta t\\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \end{bmatrix} \end{equation*}` + + Observation Model ~~~~~~~~~~~~~~~~~ +Position Estimation Kalman Filter +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + The robot can get x-y position infomation from GPS. So GPS Observation model is @@ -120,6 +201,30 @@ Its Jacobian matrix is :math:`\begin{equation*} 　= \begin{bmatrix} 1& 0 & 0 & 0\\ 0 & 1 & 0 & 0\\ \end{bmatrix} \end{equation*}` +Kalman Filter with Speed Scale Factor Correction +^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ + +The robot can get x-y position infomation from GPS. + +So GPS Observation model is + +.. math:: \textbf{z}_{t} = g(\textbf{x}_t) = H \textbf{x}_t + +where + +:math:`\begin{equation*} H = \begin{bmatrix} 1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ \end{bmatrix} \end{equation*}` + +The observation function states that + +:math:`\begin{equation*} \begin{bmatrix} x' \\ y' \end{bmatrix} = g(\textbf{x}) = \begin{bmatrix} x \\ y \end{bmatrix} \end{equation*}` + +Its Jacobian matrix is + +:math:`\begin{equation*} J_g = \begin{bmatrix} \frac{\partial x'}{\partial x} & \frac{\partial x'}{\partial y} & \frac{\partial x'}{\partial \phi} & \frac{\partial x'}{\partial v} & \frac{\partial x'}{\partial s}\\ \frac{\partial y'}{\partial x}& \frac{\partial y'}{\partial y} & \frac{\partial y'}{\partial \phi} & \frac{\partial y'}{ \partial v} & \frac{\partial y'}{ \partial s}\\ \end{bmatrix} \end{equation*}` + +:math:`\begin{equation*} 　= \begin{bmatrix} 1& 0 & 0 & 0 & 0\\ 0 & 1 & 0 & 0 & 0\\ \end{bmatrix} \end{equation*}` + + Extended Kalman Filter ~~~~~~~~~~~~~~~~~~~~~~