1. LP Problem - DMC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "026bf254",
   "metadata": {},
   "source": [
    "# Linear Programming: DMC Formulation (Simplified)\n",
    "- Author: Siang Lim, Shams Elnawawi\n",
    "- Last Updated: June 7th 2022\n",
    "- Created March 2nd 2022\n",
    "\n",
    "## References\n",
    "- O'Connor, D. L. (2007). Control Consulting Inc. DMC Training Slides.\n",
    "- Morshedi, A. M., Cutler, C. R., & Skrovanek, T. A. (1985). Optimal solution of dynamic matrix control with linear programing techniques (LDMC). In 1985 American Control Conference (pp. 199-208). IEEE.\n",
    "- Sorensen, R. C., & Cutler, C. R. (1998). LP integrates economics into dynamic matrix control. Hydrocarbon Processing, 77(9), 57-65.\n",
    "- Ranade, S. M., & Torres, E. (2009). From dynamic mysterious control to dynamic manageable control. Hydrocarbon Processing, 88(3), 77-81.\n",
    "- Godoy, J. L., Ferramosca, A., & González, A. H. (2017). Economic performance assessment and monitoring in LP-DMC type controller applications. Journal of Process Control, 57, 26-37."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "28dbbaba",
   "metadata": {},
   "source": [
    "# 1. DMC Overview\n",
    "\n",
    "DMC (Dynamic Matrix Control) belongs to a family of controllers known as 'Model Predictive Control' (MPC).\n",
    "\n",
    "> **Note:**\n",
    "One important point here is that DMC-type MPCs in the process industry are implemented with an integrated steady-state optimizer (usually LP) that generates SS targets for the MPC algorithm. This is something that's not often clear in academic papers, which usually focuses on the MPC algorithm itself. Another common misconception is that the LP is the same as an RTO layer, this is also not correct. See Godoy (2017).\n",
    "\n",
    "DMC answers 3 questions about the process and how to control it.\n",
    "\n",
    "## **Where am I?**\n",
    "\n",
    "It uses empirical **models** from step tests to **predict** how the process is moving.\n",
    "\n",
    "## **Where should I go?**\n",
    "\n",
    "It uses a steady-state LP optimizer to generate **SS targets** that keep the process at an economic optimum.\n",
    "\n",
    "## **How should I get there?**\n",
    "\n",
    "It generates a **move plan** to move the process to the optimum point based on these SS targets.\n",
    "\n",
    "# This notebook: How do we generate the steady-state targets?\n",
    "\n",
    "The LP steady-state (SS) optimizer is responsible for generating SS targets for the move calculations. We will look at this component first, and the Prediction and Move Plan modules in subsequent notebooks.\n",
    "\n",
    "<div>\n",
    "<br>\n",
    "<img src=\"./assets/img/ranade.png\" width=\"500\"/>\n",
    "</div>\n",
    "\n",
    "From: *Ranade, S. M., & Torres, E. (2009). From dynamic mysterious control to dynamic manageable control. Hydrocarbon Processing, 88(3), 77-81.*"
   ]
  },
  {
   "attachments": {
    "image-3.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAABMoAAAHkCAYAAADPbcPYAAAMamlDQ1BJQ0MgUHJvZmlsZQAASImVVwdYU8kWnluSkJDQAhGQEnoTRHqREkKLICBVsBGSQEKJMSGIWFEXFVy7iGJFV0UUXV0BWVTEXhbF3hcLKsoqFhRF5U0K6LqvfO9839z5c+bMf0pm7p0BQLuXK5HkojoA5InzpfERIcyxqWlM0hOAACqggaHAicuTSVhxcdEAykD/d3l/HVpDueKs4Prn+H8VPb5AxgMAGQ9xBl/Gy4O4GQB8PU8izQeAqNBbTc2XKPAciPWlMECIVylwlgrvVOAMFW5S2iTGsyG+BIAGlcuVZgGgdRfqmQW8LMij9RliVzFfJAZAexjEgTwhlw+xIvZheXmTFbgCYntoL4EYxgN8Mr7jzPobf8YgP5ebNYhVeSlFI1Qkk+Ryp/2fpfnfkpcrH/BhCxtVKI2MV+QPa3gzZ3KUAlMh7hJnxMQqag1xr4ivqjsAKEUoj0xS2aMmPBkb1g8wIHblc0OjIDaBOFycGxOt1mdkisI5EMPVghaK8jmJEBtCvFAgC0tQ22yWTo5X+0LrMqVsllp/hitV+lX4ui/PSWKp+d8IBRw1P6ZVJExMgZgCsXWBKDkGYi2IXWQ5CVFqm5FFQnbMgI1UHq+I3xrieIE4IkTFjxVkSsPj1falebKBfLHNQhEnRo335wsTI1X1wU7wuMr4YS7YJYGYlTTAI5CNjR7IhS8IDVPljj0TiJMS1Dy9kvyQeNVcnCLJjVPb45aC3AiF3hJiD1lBgnounpwPF6eKH8+U5MclquLEi7K5o+JU8eDLQDRgg1DABHLYMsBkkA1ErV31XfCXaiQccIEUZAEBcFZrBmakKEfE8JkAisBfEAmAbHBeiHJUAAqg/sugVvV0BpnK0QLljBzwBOI8EAVy4W+5cpZ40FsyeAw1on9458LGg/HmwqYY//f6Ae03DQtqotUa+YBHpvaAJTGMGEqMJIYTHXBjPBD3x6PhMxg2N9wH9x3I45s94QmhjfCQcI3QTrg1STRX+kOUo0E75A9X1yLj+1rgtpDTEw/BAyA7ZMYZuDFwxj2gHxYeBD17Qi1bHbeiKswfuP+WwXf/htqO7EpGyUPIwWT7H2dqOWp5DrIoav19fVSxZgzWmz048qN/9nfV58M+6kdLbCF2ADuNHcPOYk1YPWBiR7EG7AJ2WIEHV9dj5eoa8BavjCcH8oj+4Y+r9qmopMy1xrXT9bNqLF9QmK/YeOzJkmlSUZYwn8mCXwcBkyPmuQxjurm6uQGg+NaoXl9vGcpvCMI4901XPB2AkQf6+/ubvumi4dvroCnc/q++6ew/wNeEKwBninlyaYFKhyseBPiW0IY7zQiYAStgD/NxA17AHwSDMDAKxIJEkAomwioL4TqXgqlgBigGJaAMLAOrwTqwCWwFO8EesB/UgyZwDJwC58ElcA3cgaunA7wA3eA96EMQhITQEDpihJgjNogT4ob4IIFIGBKNxCOpSDqShYgROTIDmYeUISuQdcgWpBr5FTmEHEPOIm3ILeQB0om8QT6hGEpF9VFT1BYdjvqgLDQKTUQnoFnoFLQInY8uQSvQKnQ3WoceQ8+j19B29AXagwFME2NgFpgz5oOxsVgsDcvEpNgsrBQrx6qwWqwR/s9XsHasC/uIE3E6zsSd4QqOxJNwHj4Fn4UvxtfhO/E6/AR+BX+Ad+NfCTSCCcGJ4EfgEMYSsghTCSWEcsJ2wkHCSbiXOgjviUQig2hH9IZ7MZWYTZxOXEzcQNxLbCa2ER8Re0gkkhHJiRRAiiVxSfmkEtJa0m7SUdJlUgepV0NTw1zDTSNcI01DrDFXo1xjl8YRjcsaTzX6yDpkG7IfOZbMJ08jLyVvIzeSL5I7yH0UXYodJYCSSMmmFFMqKLWUk5S7lLeampqWmr6aYzRFmnM0KzT3aZ7RfKD5kapHdaSyqeOpcuoS6g5qM/UW9S2NRrOlBdPSaPm0JbRq2nHafVqvFl3LRYujxdearVWpVad1WeulNlnbRpulPVG7SLtc+4D2Re0uHbKOrQ5bh6szS6dS55DODZ0eXbruCN1Y3Tzdxbq7dM/qPtMj6dnqhenx9ebrbdU7rveIjtGt6Gw6jz6Pvo1+kt6hT9S30+foZ+uX6e/Rb9XvNtAz8DBINig0qDQ4bNDOwBi2DA4jl7GUsZ9xnfFpiOkQ1hDBkEVDaodcHvLBcKhhsKHAsNRwr+E1w09GTKMwoxyj5Ub1RveMcWNH4zHGU403Gp807hqqP9R/KG9o6dD9Q2+boCaOJvEm0022mlww6TE1M40wlZiuNT1u2mXGMAs2yzZbZXbErNOcbh5oLjJfZX7U/DnTgMli5jIrmCeY3RYmFpEWcostFq0WfZZ2lkmWcy33Wt6zolj5WGVarbJqseq2NrcebT3Dusb6tg3ZxsdGaLPG5rTNB1s72xTbBbb1ts/sDO04dkV2NXZ37Wn2QfZT7KvsrzoQHXwcchw2OFxyRB09HYWOlY4XnVAnLyeR0wantmGEYb7DxMOqht1wpjqznAuca5wfuDBcol3mutS7vBxuPTxt+PLhp4d/dfV0zXXd5npnhN6IUSPmjmgc8cbN0Y3nVul21Z3mHu4+273B/bWHk4fAY6PHTU+652jPBZ4tnl+8vL2kXrVend7W3une671v+Oj7xPks9jnjS/AN8Z3t2+T70c/LL99vv98rf2f/HP9d/s9G2o0UjNw28lGAZQA3YEtAeyAzMD1wc2B7kEUQN6gq6GGwVTA/eHvwU5YDK5u1m/UyxDVEGnIw5APbjz2T3RyKhUaEloa2humFJYWtC7sfbhmeFV4T3h3hGTE9ojmSEBkVuTzyBseUw+NUc7pHeY+aOepEFDUqIWpd1MNox2hpdONodPSo0StH342xiRHH1MeCWE7syth7cXZxU+J+H0McEzemcsyT+BHxM+JPJ9ATJiXsSnifGJK4NPFOkn2SPKklWTt5fHJ18oeU0JQVKe1jh4+dOfZ8qnGqKLUhjZSWnLY9rWdc2LjV4zrGe44vGX99gt2EwglnJxpPzJ14eJL2JO6kA+mE9JT0XemfubHcKm5PBidjfUY3j81bw3vBD+av4ncKAgQrBE8zAzJXZD7LCshamdUpDBKWC7tEbNE60evsyOxN2R9yYnN25PTnpuTuzdPIS887JNYT54hPTDabXDi5TeIkKZG0T/GbsnpKtzRKul2GyCbIGvL14aH+gtxe/pP8QUFgQWVB79TkqQcKdQvFhRemOU5bNO1pUXjRL9Px6bzpLTMsZhTPeDCTNXPLLGRWxqyW2Vaz58/umBMxZ2cxpTin+I+5rnNXzH03L2Ve43zT+XPmP/op4qeaEq0SacmNBf4LNi3EF4oWti5yX7R20ddSfum5Mtey8rLPi3mLz/084ueKn/uXZC5pXeq1dOMy4jLxsuvLg5bvXKG7omjFo5WjV9atYq4qXfVu9aTVZ8s9yjetoayRr2mviK5oWGu9dtnaz+uE665VhlTuXW+yftH6Dxv4Gy5vDN5Yu8l0U9mmT5tFm29uidhSV2VbVb6VuLVg65NtydtO/+LzS/V24+1l27/sEO9o3xm/80S1d3X1LpNdS2vQGnlN5+7xuy/tCd3TUOtcu2UvY2/ZPrBPvu/5r+m/Xt8ftb/lgM+B2t9sflt/kH6wtA6pm1bXXS+sb29IbWg7NOpQS6N/48HfXX7f0WTRVHnY4PDSI5Qj84/0Hy062tMsae46lnXsUcukljvHxx6/emLMidaTUSfPnAo/dfw06/TRMwFnms76nT10zudc/Xmv83UXPC8c/MPzj4OtXq11F70vNlzyvdTYNrLtyOWgy8euhF45dZVz9fy1mGtt15Ou37wx/kb7Tf7NZ7dyb72+XXC7786cu4S7pfd07pXfN7lf9afDn3vbvdoPPwh9cOFhwsM7j3iPXjyWPf7cMf8J7Un5U/On1c/cnjV1hndeej7ueccLyYu+rpK/dP9a/9L+5W+vgl9d6B7b3fFa+rr/zeK3Rm93vPN419IT13P/fd77vg+lvUa9Oz/6fDz9KeXT076pn0mfK744fGn8GvX1bn9ef7+EK+UqjwIYbGhmJgBvdgBASwWADu9tlHGqu6BSENX9VYnAf8Kq+6JSvACohZ3iGM9uBmAfbHbBkBs2xRE+MRig7u6DTS2yTHc3FRcV3oQIvf39b+E5htQIwBdpf3/fhv7+L9tgsLcAaJ6iuoMqhAjvDJuVHJcZhXPAD6K6n36X4489UETgAX7s/wV1VI8p6UPjeQAAAIplWElmTU0AKgAAAAgABAEaAAUAAAABAAAAPgEbAAUAAAABAAAARgEoAAMAAAABAAIAAIdpAAQAAAABAAAATgAAAAAAAACQAAAAAQAAAJAAAAABAAOShgAHAAAAEgAAAHigAgAEAAAAAQAABMqgAwAEAAAAAQAAAeQAAAAAQVNDSUkAAABTY3JlZW5zaG90EC3UewAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAAddpVFh0WE1MOmNvbS5hZG9iZS54bXAAAAAAADx4OnhtcG1ldGEgeG1sbnM6eD0iYWRvYmU6bnM6bWV0YS8iIHg6eG1wdGs9IlhNUCBDb3JlIDYuMC4wIj4KICAgPHJkZjpSREYgeG1sbnM6cmRmPSJodHRwOi8vd3d3LnczLm9yZy8xOTk5LzAyLzIyLXJkZi1zeW50YXgtbnMjIj4KICAgICAgPHJkZjpEZXNjcmlwdGlvbiByZGY6YWJvdXQ9IiIKICAgICAgICAgICAgeG1sbnM6ZXhpZj0iaHR0cDovL25zLmFkb2JlLmNvbS9leGlmLzEuMC8iPgogICAgICAgICA8ZXhpZjpQaXhlbFlEaW1lbnNpb24+NDg0PC9leGlmOlBpeGVsWURpbWVuc2lvbj4KICAgICAgICAgPGV4aWY6UGl4ZWxYRGltZW5zaW9uPjEyMjY8L2V4aWY6UGl4ZWxYRGltZW5zaW9uPgogICAgICAgICA8ZXhpZjpVc2VyQ29tbWVudD5TY3JlZW5zaG90PC9leGlmOlVzZXJDb21tZW50PgogICAgICA8L3JkZjpEZXNjcmlwdGlvbj4KICAgPC9yZGY6UkRGPgo8L3g6eG1wbWV0YT4KfnLZlgAAABxpRE9UAAAAAgAAAAAAAADyAAAAKAAAAPIAAADyAABibkY8jxIAAEAASURBVHgB7N0JvFTz/8fxz9z9Jm2WFv7ZklBEtKDVVil+SYSyZKsoslS2RJSSaFEIoZBUEmnVIjvZQoqStCmlRXfuPv/vuTXdWc6ZOTNzZubMzGseP787c9bveX5P33vPe77nexwu9RJeCCCAAAIIIIAAAggggAACCCCAAAIIpLiAg6Asxc8ADh8BBBBAAAEEEEAAAQQQQAABBBBAoEyAoIwTAQEEEEAAAQQQQAABBBBAAAEEEEAAASVAUMZpgAACCCCAAAIIIIAAAggggAACCCCAgBIgKOM0QAABBBBAAAEEEEAAAQQQQAABBBBAQAkQlHEaIIAAAggggAACCCCAAAIIIIAAAgggoAQIyjgNEEAAAQQQQAABBBBAAAEEEEAAAQQQUAIEZZwGCCCAAAIIIIAAAggggAACCCCAAAIIKAGCMk4DBBBAAAEEEEAAAQQQQAABBBBAAAEElABBGacBAggggAACCCCAAAIIIIAAAggggAACSoCgjNMAAQQQQAABBBBAAAEEEEAAAQQQQAABJUBQxmmAAAIIIIAAAggggAACCCCAAAIIIICAEiAo4zRAAAEEEEAAAQQQQAABBBBAAAEEEEBACRCUcRoggAACCCCAAAIIIIAAAggggAACCCCgBAjKOA0QQAABBBBAAAEEEEAAAQQQQAABBBBQAgRlnAYIIIAAAggggAACCCCAAAIIIIAAAggoAYIyTgMEEEAAAQQQQAABBBBAAAEEEEAAAQSUAEEZpwECCCCAAAIIIIAAAggggAACCCCAAAJKgKCM0wABBBBAAAEEEEAAAQQQQAABBBBAAAElQFDGaYAAAggggAACCCCAAAIIIIAAAggggIASICjjNEAAAQQQQAABBBBAAAEEEEAAAQQQQEAJEJRxGiCAAAIIIIAAAggggAACCCCAAAIIIKAECMo4DRBAAAEEEEAAAQQQQAABBBBAAAEEEFACBGWcBggggAACCCCAAAIIIIAAAggggAACCCgBgjJOAwQQQAABBBBAAAEEEEAAAQQQQAABBJQAQRmnAQIIIIAAAggggAACCCCAAAIIIIAAAkqAoIzTAAEEEEAAAQQQQAABBBBAAAEEEEAAASVAUMZpgAACCCCAAAIIIIAAAggggAACCCCAgBIgKOM0QAABBBBAAAEEEEAAAQQQQAABBBBAQAkQlHEaIIAAAggggAACCCCAAAIIIIAAAgggoAQIyjgNEEAAAQQQQAABBBBAAAEEEEAAAQQQUAIEZZwGCCCAAAIIIIAAAggggAACCCCAAAIIKAGCMk4DBBBAAAEEEEAAAQQQQAABBBBAAAEElABBGacBAggggAACCCCAAAIIIIAAAggggAACSoCgjNMAAQQQQAABBBBAAAEEEEAAAQQQQAABJUBQxmmAAAIIIIAAAggggAACCCCAAAIIIICAEiAo4zRAAAEEEEAAAQQQQAABBBBAAAEEEEBACRCUcRoggAACCCCAAAIIIIAAAggggAACCCCgBAjKOA0QQAABBBBAAAEEEEAAAQQQQAABBBBQAgRlnAYIIIAAAggggAACCCCAAAIIIIAAAggoAYIyTgMEEEAAAQQQQAABBBBAAAEEEEAAAQSUAEEZpwECCCCAAAIIIIAAAggggAACCCCAAAJKgKCM0wABBBBAAAEEEEAAAQQQQAABBBBAAAElQFDGaYAAAggggAACCCCAAAIIIIAAAggggIASICjjNEAAAQQQQAABBBBAAAEEEEAAAQQQQEAJEJRxGiCAAAIIIIAAAggggAACCCCAAAIIIKAECMo4DRBAAAEEEEAAAQQQQAABBBBAAAEEEFACBGWcBggggAACCCCAAAIIIIAAAggggAACCCgBgjJOAwQQQAABBBBAAAEEEEAAAQQQQAABBJQAQRmnAQIIIIAAAggggAACCCCAAAIIIIAAAkqAoIzTAAEEEEAAAQQQQAABBBBAAAEEEEAAASVAUMZpgAACCCCAAAIIIIAAAggggAACCCCAgBIgKOM0QAABBBBAAAEEEEAAAQQQQAABBBBAQAkQlHEaIIAAAggggAACCCCAAAIIIIAAAgggoAQIyjgNEEAAAQQQQCAqAlu2lsqOnS459eR0cTiisgs2igACCIQlsGePSzarNmr3bpccVSut7D/aqbAoWQkBBBBIOgGCsqSrUg4IAQQQQACB+Aps3lIqjz9dIBM+dImUimxakiO1aqbFt1DsHQEEUlagVLVDn31RJNM/KJE3Py6V7f8qihIfjnSRmtVEGv6fQ045wSHXXZUhp9XP8FmIjwgggAACqSBAUJYKtcwxIoBAUgrs2+eSio2dlhzbA13T5ImHcyzZltUb2bvXJZWaq+MsinzL1zR3yBvP5x7c0OHn5smOXQc/Wv6mT0eHjHmyfH/h7mDFd8VyVrfCcFf3Wu/yxg6ZMSnyMnlt9MCHrX+XyhOjCmTcByog87gIJSjT00q9abRZodc5bZZIJG1Wfr5LRo4rkIenqKQsP3T/VqeKDOqbKa1bZoa+MmsggAACCCSsAEFZwlYdBUcAgVQXKLvobGRNUCZZIv99niuHHGK/++Oeea5A7h7rkbpEUPGdz3HI9FfKQyJH4zyRPRFsMMiqt7V3yPOjyvcXZHHD2V+vKJbG11oTlJ12rMgP8yoY7iucGdu2l8rQZwpl9HvqYlSnqgjKwlFNvnVos0KvU9oskXDaLJfK6t94u1C6jyxWv9xCd/ddo+5RIi89linNzyUw87XhMwIIIJCMAgRlyVirHBMCCKSEQLH6+/+VKQWy7R9tnBWXrFzrkk/+UIcexrfmGtgzfdPlrt7ZtrLTbpdJb6bCrN1hFkvdSnN6bZEmJzvkmKMccm7jdGnZvPxCZ8rUQvn6h1L5eW2pfPST2ofaX9gvdWdh0zoipx7nkKNqOKRaFYc0OztdGp8V+a07O3eWypRpRbJhk0vWq/9m/KCuAveFUFJVtgsaiNSvkyYtmqZJp44qGbXgtf2fUhk+ulCenqkfkLl3QVDmlkjtn7RZJuqfNms/UgRtVkGBS67umS/vfq7aSStf6nukl+5Pl5uus9fvSSsPkW0hgAACCOwXICjjTEAAAQSSSEALlpZ/WiR9Hi+SlX+GeGCVRIo/ryDp6kLNLq9ZHxRKp3tVIhjqS13HvHxvunTtnCUVKpjrJffrmhLp1q9AVqwNdWciFzZ0yJSx2XLkEbEZh0vrLfHxJ0XS6l51P2qQELFpXZFXn86Wk060rmK14G74mEIZMV2dcCaqh6As9HMqVdagzTpQ07RZB0/5SNqsnf+65MJrnfLtuoObs/zNXZ3T5OkhOZIWm+be8vKzQQQQQACB4AIEZcGNWAIBBBBIOAHt4rN2qzzZtC20or85JEOu7mJNb6PQ9qy/dJOOefLVb/rzDKeqXGz1+9lSt07owZAWQPXo65RXF4bWE2HvN7lSsaK5QM6w3GHM+HNDqRzbVnUhVPWt97rgdJEFUytY9sTJXbtcMmJMgQybZi4gc5eJoMwtwU8jAdos2izt3IikzdLap7odnLL9H6OzzLrpD1xr33E9rTtKtoQAAgikrgBBWerWPUeOAAJJLjBzdqF07m+iu4+HQ51aIr8ttnb8Ko/Nh/T2+x+L5YwrQx+X60o1YP/bE8MfF0y7Pey0dnmy6i+TxVW9Cly/xM/sOBWIrt+qU1aVd25XT5s8/LDIuz3s3r1/QOzHp6qALIyHKhCU6dQPk/wEaLP8SExNoM1SQyOqsRFbd8mT5b+YIrNkoaUTM71u5bdko2wEAQQQQMAWAgRltqgGCoEAAghERyCjSZ6UBLk1z3fPHz2fKW1alY/j5Ts/Vp+79XLKG0tC69mllW3F1Cw5s2Fk44L99EuJNLiiwLCnlq/BN29lSaMzItun7zbNfq7VPE+2bPdfelSfdOl3e2Rj6WhPHNWeGPfYmz4BmZa9aR32TIZmBGX+9cMUfQHaLH2XYFNTvc2656F8GaXdCq7zanSCyEVN0uSkExxy4gnpcuwxaVKk2i7tISQ/rSqR9z8qlXdXqN81oX4vc4jIjkW5Uq1q7HsT6xwmkxBAAAEELBQgKLMQk00hgAACdhNoe7VT5n8XWtjU6lSRJTPi10NKM/x7W6nUaGV8S2Eg58KVFSTTgpyvzwCnjHvPnN3Aq9Nk2CM5gYoVlXna7Wrp9dXDDnyvD1WHOqd6imlOTmQXcNrDDroPPtArsaLI7W0ccnn7DGlxXqZkqFxwzvwi6XC3uuLUedKl5wETlHlq8D6QAG1WIJ3A81K1zfrxp2I5vYtKuXya65OOVk++fMrclxj5+S55bmKh3PuyasxCeCDOnZenybNDY9/2Bz4TmIsAAgggEKkAQVmkgqyPAAII2Fjglruc8tI8n6sHE+X9cUaWNDg1Pj2ktOINGpovQ173TX9MFFxdr7i+tybk0wasP6ylumIy02tKhUglX1WI+eDOP6wslobaBaLP64lb0uSBeyK/eJu7sEjGvFIsva9Ll0suztI9vut6O2Xy4sDnGEGZTwXx0VCANsuQJuiMVG2zzr0sTz5b7c1zW3uHjB2eG/KXJnv2uOTym53y0Y/e2zP8pDrtFnydK1lZkX0pYbh9ZiCAAAIIxEWAoCwu7OwUAQQQiI1AKD0MPEt0dUuHvPlC+ON8eW4r1PeFhS7JbuJUXaJCXVMtrwIr1zfWBGXa3gc/mS+PvmousHv/2Uzp0NaCrmwhHPYDj+XLMO22SM+XuiVyx/IcqVYt8rHJPDdr9P6tdwrlmocDj4VHUGakx3RfAdosX5HQPqdam7VkWZG0uc3724wbLnTIpLHh//7Seupq5+H49wN/AeCumbcez5CuV9jnITjucvETAQQQQCB8AYKy8O1YEwEEELC9wL0P58vT7/gEKWZKrb4c/2tRjhx9VGzCFs8ivfZGgdwwxPtevvTKqseW1nEqSHiWVUV9u/+FdUGZ1rug8rlqp97XYZ7FPfhee1rbwret2/fBDQd4U/mcPNmz03uB6y9wyKvjwr9I9N5a8E+LlxbJ+T0DAxGUBXdkif0CtFmRnQmp1mZdqXp/vfNJeaCljUf21fuR9+7VwrKLuuaZ6lnW/BSRj2fGtu2P7CxhbQQQQACBYAIEZcGEmI8AAggksMB9g/Jl5LQwgjJ1zHd1TpNnnoj89r1Q+Y5vnSd/bPFea/Rd6XLniyo8U8NxBXrlVlWLfG7tBUsot4H+9mG21DleG+U++q+ycXmu8L/tMtYPFvjqm2Jp0s2/HJ4CBGWeGrwPJECbFUjH3LxUabO0p/FWOUd9keH+XkWNFrBuTrYcd4w1bXBZ6NhSbT/IFzSivpdwfWft7x1zNc1SCCCAAALREiAoi5Ys20UAAQRsIKB30XnasSI/rjdROHUX4e5Pc6VSpdiNvbLskyJpdbNP7yRVjv++yJWKzdXVShyCsrJxf5qrscrcF2MB6Hp3dMhzT8WmN9ddD+TL6JneIehxNdWF4pLYXrARlAU4IZgVsgBtVshkfiukSpv14qsFctuT5Q3ziF5pct+d1n65M0o99feeceX78MM+MMH5feQPTzHaNtMRQAABBGIvQFAWe3P2iAACCMRMQO+ic8aIDOncP/CYUu4CPnlbmgzoZ+2Fh3vbej8vudYpH64ov41GW+beK9PkqcdyxHGmSsniEJRpZbi1n1MmzvUulzbd76V6NOxS4WLlytENFwsKXJLTWAWHBd4lGHdPutx+ixpdOoYvgrIYYqfArmizrKnkVGizOnZ3ygdfH2iX1fiU+artzc62tu3VnoaZq7W1gTvNyi/vZcvJJ1nTk82aM4CtIIAAAghEIkBQFoke6yKAAAI2F9C76Cz6qYIc2TJP/t1hovCHqOsDNeZXZgzGqF+3vkROaKeSH888Sl3zbPwoR46qlRbXoGzN7yVyUgefVMqAb5gKFwdGOVx8c1qhXDvIJ+xU12i7P49tD0CNgKDM4ERgclgCtFlhsfmtlOxtVonq5JXRUH1zcqAD8ksPpMtN10XnS4Irejhlxmeev5j8uGXuuExpe0EMflH675opCCCAAAJRECAoiwIqm0QAAQTsIqB30Vn8cwV5ebL3LSuByjvpoXS5oVt0LkA899t3oFPGzvK+GLmyuUPenrj/VsZ49ijTytnuaqfM+867fJ7lP/hehYtFX1aQDNW7LFqvxh3y5OvfvbferbVDJk+IzW2fnnsmKPPU4H2kArRZkQqWr5/MbdYXXxVLs+sOdPNSd5sXqjY3Wl/ovDqlQG58PPDtl/PHZ8pFbQjKys8+3iGAAAKJLUBQltj1R+kRQACBgAJGF53FxerWvabqdhL1v2CvmkeIbF4e3XGv9u51SSVtUGaf4cm+fTtLzjh9f+IU76BMd/w0A7zJgzOkW9csg7mRTf51TYmcfKl/77bPXsuSZk2imM4ZFJugzACGyWEJ0GaFxaa7UjK3WUOfzpcHJ+4fo/HRHmkyqH/0hgj49vtiadQ18L2Xa9RDBE48gVsvdU9EJiKAAAIJKEBQloCVRpERQAABswJGF53p6u/5IU/ly6CXvQeDN9runDGZ0v6i6H1b/uz4Auk3xvsb+yZ1Rb6YXR7QxTso02xOOj9P1mwyUiqfXqeWyG+Ly8tePifyd3feny9j3vWut+oqzNwa5TDTqOQEZUYyTA9HgDYrHDXjdZK1zVr9W4m8816RrPjJJROfzpbDD0szRohwzsZNpfJ/56sHugR4Fa6MXo+2ALtlFgIIIIBAlAQIyqIEy2YRQAABOwgEuujc+a9LDtOeJOkz1JVeuX1DK71lwp1WqjKfrGZ5UrLbewuzn8mQju3Ke2XZISjTHRvMu9gHPy15MVNatbA2XDQaxH9Un3Tpd3v0b489eHAebwjKPDB4G7EAbVbEhF4boM3y4gjrw9/bSqVGiwBBmbrd3rUiOl+MhFVgVkIAAQQQiFiAoCxiQjaAAAII2Fcg0EWnVurb73PK+PdNjLullv36zSw560zrb+17b06h/O8e77SuUjXt6ZEVxOHxADM7BGXFqpiZTdUA0v8Fr/MLGzpkwVRrxwzTvehVHSl2qKe9VavqgRW8eJYtQVBmGSUbUgK0WdaeBrRZkXsGa+POOE7k27kEZZFLswUEEEDAPgIEZfapC0qCAAIIWC4Q7KLzjz9L5Pi2Pk+aNChFp2YOmTnJ2uBH21XTS/PkyzXeOx1/X7r0usm7h5QdgjKtlE+Nzpf+E7xvffQuffmn1R9kS9061o1bozeIf5fzHDLtJevrpfwoAr8LdhGprb1pSY7Uqhm9W6MCl5C5iSRAm2V9bdFmRWaq+wWFxyZvvMghr4yJXxvsURTeIoAAAghYJEBQZhEkm0EAAQTsKBDsolMrc+cbnTLzcxO9ylSHpXXzsuW4Y6wLfn5YWSwNu/gMkqzysbwvciU317uHlF2Csj17XFJZe/CAdyc43eq/pZ1DXnzGmgsoo0H8l07MlJbNrb3FU/dgDCYSlBnAMDksAdqssNgCrkSbFZAn6MzBT+bLo68afzmy9CXVBp8XvzY46AGwAAIIIIBAyAIEZSGTsQICCCCQOAJmLjpXfFcsZ13tE1YZHGKvDg4ZP9Ka4EfbRfdeTpmyxDuku/+aNBk6yP8JZnYJyrRy9x3olLGzvMutTfd7qUxxx/IcqVYt8t5UeoP4V6wqsvfz+N7yQ1DmV+tMiECANisCvACr0mYFwAky69qeTnlzqX57b4c2OEjxmY0AAgggEIYAQVkYaKyCAAIIJIqAmYtO7Vj0bn/UPUYLg59t20uleks1QLLnF/UqT9qibtOrUd0/WLJTULZ+Q4kcd7G5W1Yf7ZEmg/r7B3+6vgYTjQbxH94zTfrfFdm2DXZpejJBmWkqFjQhQJtlAimMRWizwkA7sEqgJ4cOuy1NBvaLbxsc/pGxJgIIIICAkQBBmZEM0xFAAIEkEDB70fnhgiK5pG+RqSMefGOaPDIg8guDR4bly2OveaZkIt1aO2TyBP0ea3YKyjSoTjc4ZdYX+r0MvCBVh68CdStpVpb3raReywT58MbbhdLtEZ97PVWW+PeyHDnyCP9QMcjmLJ1NUGYpZ8pvjDYreqcAbVbotmt+L5GTOqgvRfReqknfslT/ix29xZmGAAIIIJA4AgRliVNXlBQBBBAIWcDsRae24VrN82TLdhO7UBmZUwU/OTnhBz+FhS7JbqLG+VL/83z9OCNLGpyq/2RNuwVlX35dLE27m7tlddJD6XJDN++HE3ged7D3eoP4X9bEIbNe0w8Vg23PyvkEZVZqsi3arOidA7RZodv2fyRfnnrb+wsd91baN3LInDfi3wa7y8NPBBBAAAHrBAjKrLNkSwgggIDtBEK56Jw0pUB6PF5i6hgm9E+Xnj3CD35ef7NQrn/Mu4fUeSeLLH/XeLwtuwVlGlTDdnnywx/ByWrXEPlzqfGxBdqC0SD+88dnykVt4j+ANEFZoNpjXqgCtFmhioW2PG2Wea8S9esw4+w89XQZnXVUR95Vs7KlXl3rHm6jsxcmIYAAAgjESYCgLE7w7BYBBBCIhUAoF51F6s7LrCYGFwU+hdUGMN7zWQVxhNmprE6bPFm72Xujc8dlStsLjIMfOwZlM2cXSuf+3oGf91GVf1o4IVMuaG18fOVLer/TG8RfKqmh3b4M3997D5F9IiiLzI+1vQVos7w9rP5Em2VedPaHhXLZ3frt+yM3pMnggZEPQWC+NCyJAAIIIBBLAYKyWGqzLwQQQCDGAqFcdGpFe/KZfLn/Bf3bTHyLPvOpDOnUMct3ctDPyz8tkhY3eY+HVv0Ika3LA/e4smNQVqqosprlScnuoIctreuLLJ4e+Bh9t2I0iP9jN6XJw/fZ4yKNoMy31vgciQBtViR6wdelzQpu5F6i/TVOmfut/ziU2u+rvxZXkMzQv/dwb5qfCCCAAAI2FyAos3kFUTwEEEAgEoFQLzp373ZJlXPVwGH6X6J7FeW0Y0V+mBda8KNtoEM3p8z5xvviY+L96XLz9YFv5bRjUKYdz5jnC+TOZ83dsvqzulXnlHrmb9XRHcRf9eLbrJ4MWrNGfAfx145dexGU7Xfg/60RoM2yxjHQVlK9zQpk45639o8SqdNe/8nGn7+eJU0b64+l6V6fnwgggAACiS1AUJbY9UfpEUAAgYACoV50ahvTvdXPYC+fvpol5zQ1f8Hwx58lcnxbn4sP1TEq/8tcyc4OfB+nXYOyfftcUrGZChdNjOt//QUOeXWc+cGf9Qbxb3uGQ+a+ZX4bBlVn2WSCMsso2ZASoM2K/mmQ6m2WGeE2V+TJkp/8l+zVwSHjR9qn/fUvIVMQQAABBKwQICizQpFtIIAAAjYVCOei86+NpVL7wnwR705fukfY7kyHfPim+YsGvRBu0PVp8uj9wW8jtGtQpsEMGJwvI6aauGVVdSb7e2mOHHlE8N5gRoP4fzA6Uy652D73/BCU6f7TYGKYArRZYcKFuFoqt1nBqGa8VyhXDPDvVp2rxubctjBXKlYM/KVOsO0zHwEEEEDA/gIEZfavI0qIAAIIhC0QzkWntrOrbnHKtOUmkjK17K/vZ8tJJwa/nfC//1xyqNbzynN4MrXatmU5csThwYMjOwdlmzaXytEXqHDRRFb28HVp8tgDwYNBvVBRKoqUfFVB0oJzhX3OhLoiQVmoYiwfSIA2K5COdfNSuc0KpFjW266F+j21z2cp1XH6lxnZcvJJwX/X+azJRwQQQACBBBQgKEvASqPICCCAgFmBcC86f1hZLA27mLiXUBXkxosc8sqY4L3KRk8okLtGe4/lZXZd7XjtHJRp5bvmNqe8tcxEuGjiVlOjQfwf6p4mQx4MHrJp5YnVi6AsVtKpsR/arNjVc6q2WYGE73koX0ZN9//GY9bTGXLZJaE/vCbQvpiHAAIIIGBfAYIy+9YNJUMAAQQiFgj3olPbcfNOefLJKhNFUL2btqrbCasfadzNSXvSWu45eVK4y3t7v7xn/ht6uwdl3/9YLGdcaS5cfGFgutx6g/HDC3QH8Vd0fy7KkdpHGzt768bmE0FZbJxTZS+0WbGr6VRts4yEf15VIvUv9xlDUy089NY0uf9ue31BYXQMTEcAAQQQsEaAoMwaR7aCAAII2FIgkovOBYuL5OLenvdJGh/i/dekydBBxhcSsz8slMvu9h7zpU0DkY/eMf/UTLsHZZrOuZflyWerjZ3cc6ofocLF5cbHrjeIf6he7n1F+ydBWbSFU2v7tFmxre9UbLP0hAsLXXJmB6f8vMF77u2XOmTciOA9pr3X4hMCCCCAQKILEJQleg1SfgQQQCCAQCQXndpmj2mVJxu2BtiBe5a6I2XfF7lSoYL+IMfNLs2TL9a4F97/86PnM6VNK/OD0idCUPbhgiK5pK+5cPHDsZnS7kL/4zcaxH/mUxnSqaP9bv0hKPM+r/kUmQBtVmR+oa6dim2WnlH3Xk6ZssT71vnO5zjknZdzxaH/a01vM0xDAAEEEEgSAYKyJKlIDgMBBBDQE4j0onPK1ELpPti7J5jefrRpz96ZLnf28r+d8MefiuX0K7xvSaxdQ91GuNS4R5XePhIhKHOp66wjzsuTHTv0jsB72nkniyx/199AdxB/1aGh6OsKkqEGlLbbi6DMbjWS2OWhzYpt/aVim+UrPHJMgdw33nv8zFaniiyYWkEy/b/L8F2dzwgggAACSShAUJaElcohIYAAAm6BSC86i1VGltk0T+Q/9xYD/KwkUvx5BUn3eSjYdb2dMnmx9zf1rz6cLtdf6x+qBdi67Qfzd5f9xVcL5LYnvS+63PN8f/44I0sanFqefpUN4n+2euKad64o/bumyfDBxre2+m43lp8JymKpnfz7os2KfR2nWpvlKTxnfpF0uEv1Avb4FdX4RJFl7+RKTg5dyTyteI8AAgikkgBBWSrVNseKAAIpJxDpRacGpvdtuxHkW49nSNcrym8P3La9VKq3zBfxfIiY6kRV+GXo39QnQo8yzSU/3yW5TVXYpQ472Ova1g6ZMqF8/BujQfzXzsuW44/1SSCDbTxG8wnKYgSdIruhzYp9Radam+UW/uXXEjm1ixq83+Nu+VNqi3w5M1cqViQkczvxEwEEEEhFAYKyVKx1jhkBBFJGwIqLzr17XVLpHBX8eFxMGAHWPUpk9UfltxMOfjJfHn3VMyUTGXJzmjx0b+i9oxIlKNNsBg3NlyGvex+3rpl6gOWWJTlSo/r+J1medUmerFjrveS59UQ+mVVu6j03/p8IyuJfB8lUAtqs+NRmKrVZmvDOnaVyZLt8Kdld7n2sGhLgu1m5UqUKIVm5Cu8QQACB1BQgKEvNeueoEUAgRQSsuOjUqPS2Y0S4+IVMad0yU7SniGU3UQGb+t/Bl+oUteOTXKlWNfQLkUQKysp60rVSXcpM3IHpfmLoqtUlcsplqneDz2vqExlyVefyXno+s+P+kaAs7lWQVAXQa2uKf/a/pTvYQettx2gd2iyRVGqztB5053V2en0pUfUwkV9n5ciRR+z/0sLoXGE6AggggEBqCBCUpUY9c5QIIJCiAnoXi+FcdG7aXCpHX+BzC6WBaev6IounV5DJbxXKdY96Pwjg1vYOeWFU+a2GBpvQnZxIQZl2ADf2ccqrCz0GvtE9KjVRDdWWp54YOuDRfBk7y2d5Na/wm9BvUzXaVTSmE5RFQzV1t0mbFb+6T4U2Sxt3s921ebLoBw/nyiIb382Ro2oRknmo8BYBBBBIaQGCspSufg4eAQSSXcCqi07NqVsvp7yxxCfIMQDUBqm/vE+h/L7Ze4E1c7LlxBPCG2sr0YIyox5i3iL7Pz19R7rc86LqfuYziP9dndPkmSdCv01Vbx/RmkZQFi3Z1NwubVb86j3Z26xSdTd8l5ucMvNzj99jh4ism5ktxx0T3u+l+NUWe0YAAQQQiKYAQVk0ddk2AgggEGcBKy86f/qlRBpc7n9roN4hpqtv6D3HftGWaXuGQ+a+FV5vMm39RAvKtDJfcGWefPSj9i6816/vZ8tJJ9r7Ao6gLLy6ZS19AdosfZdYTU3mNqtHX6dMWuARkqnvIH59x/5tbKzqnv0ggAACCJQLEJSVW/AOAQQQSDoBKy86NZzzu+TJ4pXhMS19KVNanpcZ3spqrUQMyhYvLZLze5p4CoKOSqMTRL6ZY99B/N1FJihzS/DTCgHaLCsUw99GsrZZ/R7Ml2dneDxgRQ37+MObWXJa/YzwsVgTAQQQQCBpBQjKkrZqOTAEEEBA5J6H8mXUdI+LA4VSuDL8Ma+WLCuSNreFHvz4Pg0znLpJxKBMO87jWuXJ+q2hH/Hrj2RI96vtO4i/+4gIytwS/LRCgDbLCsXItpFsbdajw/Nl8CSP34MqG/vqtSw5u5H1IZl2e+fsDwullfpSiKdnRnYesjYCCCAQTwGCsnjqs28EEEAgygJ39HfKc7M9bjVR+/tvRa4cckjoT510F/XENnl+Y4+55xn9fOOxDLnmyshCn0QNyvQeamDkdHC6onJ+lSs5OeHX08FtRfkNQVmUgVNs87RZ8a/wZGqznh1fIP3GeDx+WN3J/vGLmdL83PB7NweqoRHP5suA50tFe4rmzk/t3yM40LEwDwEEEEhlAYKyVK59jh0BBJJe4La7nfLih95B2Y7PcqRatfCf7jV1eqFc/ZD30ywDQlYUKfqigmRE+OV9ogZlRaoDXlaTPPVoy4BKXjN7d3TIc0+FP56b18ai/IGgLMrAKbZ52qz4V3iytFmTphRIj8c9QjL1a2/e2Ey5+PzohGRl43heocbxVL3K7rjMIWOHJ0YbHv8zjhIggAAC9hMgKLNfnVAiBBBAwDIBvSdVrl+QI8fUDj8oK1HXHRnNVOqzx1wxh/dMk/53Rf7kRkdDtc/8wPvMqiJSoEI5u72GjcqXB170uPUnSAFXqqew1T/F3oP4uw+BoMwtwU8rBGizrFCMfBuJ3mZNn1UoXe5XX+i4vydSnXPffSpD/tchsp7NRrL//eeS49s6Zfs/+5f45b1sOfmkxGjDjY6J6QgggEAqCxCUpXLtc+wIIJD0Aq0758nSn70P8zM1NkuzJpF17/K7ncV7F+Wf1G52fZorlStHfguho4EKyoINj6Z6r7m+sV9QtmuXS6qe5xQx0RHv1NoiPy2w3zGUV6r3O4Iybw8+RSZAmxWZn1VrJ3KbNW9RkbTrq35ZeHw3YcXt/0a2e/e65LwuTvlx/f4lzjxeZMWHidOGGx0X0xFAAIFUFiAoS+Xa59gRQCDpBf6vZZ5s/Nv7MKc+kSFXdY7sW/V9+1xSsZkKfgq9t+376fZLHTJuROS3n2gDJKefqoIyd+8A3x25P6uOa67v7XmB0vtep0z4INgBiLz8YLr06J7tPiLb/1z2SZG0ujlwgrlpSY7Uqhl+L0bbI1BAywRosyyjjHhDidhmffJZkTS/RbVHHndcTrw/XW6+Pjpt6rffF0u3+wpl1V/l3JMHZ0i3rpH9ji3fGu8QQAABBOIhQFAWD3X2iQACCMRAID/fJblnqjDL41t1bbd3dU6TZ56I/FbI+x/Nlyff8tm453GpTmTr5mXLccdEfvvJps2lcnSbIPddHth3JE/19Cy+1e/X/lEiddqr8WsCZWWqB94+NYh/hQqR98CzuvxG23tvTqH8757AXeU2Ls6Ro2oRlBkZMn2/AG2Wvc6ERGuztNCqUXf17Y1Hbl+pmsjj10f+O8hdM9qXNv/udsm2f1zyyUqXrPzTPefAT5WP5X+dK9nZidOG+xwBHxFAAAEElABBGacBAgggkKQCi5YUyYW9PK4YDhznsTVE/lgaea+rrX+XSs3WKrwyyMo6NnbI7Ncj702mFXvuwiJp38f/WPSq7tu3s+SM0yO7tVRvu1ZM69DNKXO+MU7Kbm7rkInPWmNmRXnNbOPpsQVy73Me3Td0VloxNUvObGjPOtEpLpPiJECbFSf4ALtNlDZr1eoSOeVK9UWE+l88X33+55AxTyZWGx5PL/aNAAII2FWAoMyuNUO5EEAAgQgF2lyRJ0t+0t/Ioucz5fxWkT/564Y7nPLaIv3gx4qx0Nyl13sSnnue78+rWzrkzRfseaHy6efFct6NxverJmKg1O5qp8z7Tv8ccNfN9OEZ0vkybkVye/BTX4A2S98lnlMToc1av6FEjuukErJ98ZTav+9Vs7OlXl3rerDF/4goAQIIIJCaAgRlqVnvHDUCCCS5wPiXCuT2kca9fCpWFfltVo7UqB7Z7XC/rimRky/1/wq/wTEiP86PvNeaVk0rfy6W07qocMmg55pfVao7Xn7/MFtOOM6eFyv1L8qTnzf4lVrq1FJ1stgaM/+tR2fKmt9L5KSOqv4D52TSq4NDxo+0Z3gZHRm2GqoAbRZtVqjnjLb8lq2l8n+X5UvJ7nDWtnads+uIfPVBYrXh1gqwNQQQQCB5BAjKkqcuORIEEEBAfltbIo+PKpTXPwqSXGhWlUWmDcyQK/6XJY4IhlNpq3oUzffpUTRtWIZ06RRZDyJtvKK3phdKj+Eq8DN312X5GaCuVaY/kiGXXZIlGTa742/azEK56gH/Mb3G35cuvW6KzoDT5TDWvftrY6k0uyZfNm0zt80ZIzLk8ksjOyfM7YmlEkmANutAbdFmhXza7txZKvVVSLZle8irRmWFaD5ZMyoFZqMIIIAAAoYCBGWGNMxAAAEE7C2gXSS8+0GRbNjkkt/Wu2TWjy5x/htGmVVHnw71HXLmKQ6pVcMhbVpkyIknmO+NVfaUsR4eSVYlkeLPK0i6+U1IicrC3nqnULZuc8kGNXD/ilUu+WydOhaPzYZxZCIqd7qgnsiZJ6fJ0TUdcvYZ6dK0cXyTM+1YM5qpJ3ju8TgiZbXni1w59NAIEkuPzVn91qVyV+1JpxtV3Wi9COcsKpGXFqqJxp0W/YugDu0yNW7dpRekyRmnpUvNGmlyxOFpIZ0n/htlSiIJ0GaZqC3aLBNIIv/955LGnZxeT5s0tWK0FlL1VqAG8c/KsmcbHq3DZrsIIIBAsgoQlCVrzXJcCCCQ9AKffVEs595gPN5VuAAjb0+Xe/qE1rPpVHU74S8Hbid8+o50ufuO0NbXLnoOPUs9oTPKr87nOGT6K/G/BXDUuAK5Z1x5ytS9jUNeHx//cunxX3hVniz6QW+ORdPUdeUt6iEGLz5jz+O36CjZjBKgzQr9NKDN0je7o79Tnpttoue0/uqWT73z8jR5dmjkT5O2vGBsEAEEEEAgLAGCsrDYWAkBBBCIv8Dv60rkxPb+44NFWrJwgrLpswqly0B1O6F6PsDez3OlYsXQvlWPVVB2x2UOGTs8/oFM2fE2U8HggR5zX0zOkiZnx7enm9F5c9UtTpm2PLoXpHYOCo1cmB66AG1W6Ga0WfpmAwbny4ipZgeu1N+GlVNXf5AtdeuE0I3ayp2zLQQQQAABywUIyiwnZYMIIIAAAggEF7j34Xx5enqp1K4u8udSBoAOLsYSCCAQTwHarHjqs28EEEAAgVgKEJTFUpt9IYAAAggggAACCCCAAAIIIIAAAgjYVoCgzLZVQ8EQQAABBBBAAAEEEEAAAQQQQAABBGIpQFAWS232hQACCCCAAAIIIIAAAggggAACCCBgWwGCMttWDQVDAAEEEEAAAQQQQAABBBBAAAEEEIilAEFZLLXZFwIIIIAAAggggAACCCCAAAIIIICAbQUIymxbNRQMAQQQQAABBBBAAAEEEEAAAQQQQCCWAgRlsdRmXwgggAACCCCAAAIIIIAAAggggAACthUgKLNt1VAwBBBAAAEEEEAAAQQQQAABBBBAAIFYChCUxVKbfSGAAAIIIIAAAggggAACCCCAAAII2FaAoMy2VUPBEEAAAQQQQAABBBBAAAEEEEAAAQRiKUBQFktt9oUAAggggAACCCCAAAIIIIAAAgggYFsBgjLbVg0FQwABBBBAAAEEEEAAAQQQQAABBBCIpQBBWSy12RcCCCCAAAIIIIAAAggggAACCCCAgG0FCMpsWzUUDAEEEEAAAQQQQAABBBBAAAEEEEAglgIEZbHUZl8IIIAAAggggAACCCCAAAIIIIAAArYVICizbdVQMAQQQAABBBBAAAEEEEAAAQQQQACBWAoQlMVSm30hgAACCCCAAAIIIIAAAggggAACCNhWgKDMtlVDwRBAAAEEEEAAAQQQQAABBBBAAAEEYilAUBZLbfaFAAIIIIAAAggggAACCCCAAAIIIGBbAYIy21YNBUMAAQQQQAABBBBAAAEEEEAAAQQQiKUAQVkstdkXAggggAACCCCAAAIIIIAAAggggIBtBQjKbFs1FAwBBBBAAAEEEEAAAQQQQAABBBBAIJYCBGWx1GZfCCCAAAIIIIAAAggggAACCCCAAAK2FSAos23VUDAEEEAAAQQQQAABBBBAAAEEEEAAgVgKEJTFUpt9IYAAAggggAACCCCAAAIIIIAAAgjYVoCgzLZVQ8EQQAABBBBAAAEEEEAAAQQQQAABBGIpQFAWS232hQACCCCAAAIIIIAAAggggAACCCBgWwGCMttWDQVDAAEEEEAAAQQQQAABBBBAAAEEEIilAEFZLLXZFwIIIIAAAggggAACCCCAAAIIIICAbQUIymxbNRQMAQQQQAABBBBAAAEEEEAAAQQQQCCWAgRlsdRmXwgggAACCCCAAAIIIIAAAggggAACthUgKLNt1VAwBBBAAAEEEEAAAQQQQAABBBBAAIFYChCUxVKbfSGAAAIIIIAAAggggAACCCCAAAII2FaAoMy2VUPBEEAAAQQQQAABBBBAAAEEEEAAAQRiKUBQFktt9oUAAggggAACCCCAAAIIIIAAAgggYFsBgjLbVg0FQwABBBBAAAEEEEAAAQQQQAABBBCIpQBBWSy12RcCCCCAAAIIIIAAAggggAACCCCAgG0FCMpsWzUUDAEEEkFg165dUlRUFNOiZmZmSpUqVWK6T3aGAAIIIIAAAqEL7N69WwoLC0NfMcI1KlWqJNnZ2RFuhdURQACB1BQgKEvNeueoEUDAAoGlS5dKmzZtxOVyWbA185s455xz5NNPPzW/AksigAACCCCAQFwELrzwQlm0aFHM9z1t2jTp0qVLzPfLDhFAAIFkECAoS4Za5BgQQCAuAh9++KFccsklMd83QVnMydkhAggggAACYQkQlIXFxkoIIIBAXAUIyuLKz84RQCCRBQjKErn2KDsCCCCAAALRFyAoi74xe0AAAQSsFiAos1qU7SGAQMoIEJSlTFVzoAgggAACCIQlQFAWFhsrIYAAAnEVICiLKz87RwCBRBYgKEvk2qPsCCCAAAIIRF+AoCz6xuwBAQQQsFqAoMxqUbaHAAIpI0BQljJVzYEigAACCCAQlgBBWVhsrIQAAgjEVYCgLK787BwBBBJZIFBQduqpp0rdunWjcnj16tWToUOHRmXbbBQBBBBAAAEErBN49913ZfLkydZt0GNLX3zxhWzZssVjSvlbnnpZbsE7BBBAIFQBgrJQxVgeAQQQOCAQKCgbNmyYDBw4ECsEEEAAAQQQQCAqAmPHjpW+ffvqbpugTJeFiQgggIApAYIyU0wshAACCPgLEJT5mzAFAQQQQAABBGIjQFAWG2f2ggACqSdAUJZ6dc4RI4CARQIEZRZBshkEEEAAAQQQCFmAoCxkMlZAAAEETAkQlJliYiEEEEDAX4CgzN+EKQgggAACCCAQGwGCstg4sxcEEEg9AYKy1KtzjhgBBCwSICizCJLNIIAAAggggEDIAgRlIZOxAgIIIGBKgKDMFBMLIYAAAv4CBGX+JkxBAAEEEEAAgdgIEJTFxpm9IIBA6gkQlKVenXPECCBgkQBBmUWQbAYBBBBAAAEEQhYgKAuZjBUQQAABUwIEZaaYWAgBBBDwFyAo8zdhCgIIIIAAAgjERoCgLDbO7AUBBFJPgKAs9eqcI0YAAYsECMosgjywmc2bN8vKlStF+7l7927ZtWtX2c89e/ZIZmamVK1aVapVqyaHH364nH766dKgQYOy6VaWYu/evTJkyBBxuVwyaNAgOfTQQ63cPNtCAIEkEqDNSqLKTNBDIShL0Iqj2AggYHsBgjLbVxEFRAABuwoQlIVfM//9958sWrRIli5dKj/++GPZfzt27Ahpgzk5OXLmmWfK5ZdfLtdcc43UrFkzpPV9F9YCuYsuuki+/PLLslnjx4+XXr16+S7GZwQQSEEB2qwUrPQEOGSCsgSoJIqIAAIJKUBQlpDVRqERQMAOAgRlodXC+vXrZfbs2TJnzhxZtmyZFBQUhLaBAEunpaVJ+/bt5dFHHy0LzwIsqjtLK0vbtm3Lgjv3AlOnTpWrrrrK/ZGfCCCQYgK0WSlW4Ql4uARlCVhpFBkBBBJCgKAsIaqJQiKAgB0FCMqC10pxcbG8//778vzzz8vChQvLbmkMvlZkS3Tp0kW03mDaLZpmXqWlpXLllVfKjBkzvBbXbquKtJea1wb5gAACthegzbJ9FVFADwGCMg8M3iKAAAIWChCUWYjJphBAILUECMqM6/uff/6RcePGycSJE8vGHDNe0n+Ow+GQevXqybHHHiuHHXZY2X+VKlUSbfywLVu2lP23bt062bhxo//KB6YcffTRMm3aNGnWrJnhMu4Z2u2VWpDn+apbt66sXr3acxLvEUAgiQVos5K4cpP40AjKkrhyOTQEEIirAEFZXPnZOQIIJLIAQZl/7WnjfD399NPyzDPPlAVb/kv4T9Fum2zZsqW0adNGmjRpIo0bN5bKlSv7L+gz5fvvvy/rBTZz5kz55ZdffOZK2UD/Tz31lNx5551+89wTBg8eXHa7pvuz++ctt9wiL774ovsjPxFAIEkFaLOStGJT5LAIylKkojlMBBCIvYB6shcvBBBAAIEwBNRYWy7Vauv+N2zYsDC2mLir5OXluUaMGOFSPcB0PfSczj77bJcK1FzqFseID1yrixNOOEF3371799bdvro9U3d5raxTpkzRXYeJCCCQHAK0WclRj6l+FGPGjDH8PaZ6Vac6D8ePAAIIhC1Aj7LYZ5PsEQEEkkSAHmX7K/Kjjz6SW2+9VbTbIYO9tN5jnTt3lgceeEAaNmwYbPGQ5ufn58vw4cNFhZR+Dwp44oknyvbp3uD06dPLBurXxifTe/3111+i3b7JCwEEkk+ANiv56jRVj4geZala8xw3AghEW4CgLNrCbB8BBJJWINWDsn///VfuvfdeeeWVV4LWsRaQdevWTe6///6y8ceCrhDBAgsWLJD//e9/4nQ6vbby+uuvS/fu3WXJkiXSrl07vzDNvbDqmSa///67+yM/EUAgSQRos5KkIjmMgwIEZQcpeIMAAghYKkBQZiknG0MAgVQSSOWgTN3qKDfffLNs3bo1aJUff/zx8uqrr0rz5s2DLmvVAosXL5aOHTuKur3q4CZzc3Pl5Zdflp49e4o2LpHRq0ePHmXLGc1nOgIIJJ4AbVbi1RklDi5AUBbciCUQQACBcAQIysJRYx0EEEBACaRqUDZy5EgZMGCAGN226HlyaLdkaoP7V6xY0XNyTN5rPcvatm0ranCCkPb32muvyXXXXRfSOiyMAAL2FaDNsm/dULLIBAjKIvNjbQQQQMBIgKDMSIbpCCCAQBCBVAvKioqKpFevXqZ6W9WqVatsOS2oiudLDeQvEyZMCKkI69evl2OOOSakdVgYAQTsJ0CbZb86oUTWChCUWevJ1hBAAAG3AEGZW4KfCCCAQIgCqRSUaWP7aIPwa+N7BXtdffXV8txzz0nVqlWDLRr1+Xv37pVTTz1VtMH5zby0gEwLynghgEBiC9BmJXb9UXpzAgRl5pxYCgEEEAhVgKAsVDGWRwABBA4IpEpQ9ttvv0mHDh1kzZo1Aes+OztbJk2aJFpQZqfXrFmzpFOnTqaKdP3115eNp2ZqYRZCAAFbCtBm2bJaKFQUBAjKooDKJhFAAAElQFDGaYAAAgiEKZAKQdnSpUvLepLt3LkzoNKhhx4qs2fPllatWgVcLh4ztbHU6tSpI3/88UfQ3WtP8LzxxhuDLscCCCBgTwHaLHvWC6WKjgBBWXRc2SoCCCBAUMY5gAACCIQpkOxBmRYaaU+I1Mb5CfQ64ogjZO7cudKoUaNAi8V13tChQ+XBBx8MWoa1a9eK9pROXgggkHgCtFmJV2eUODIBgrLI/FgbAQQQMBIgKDOSYToCCCAQRCCZgzLtFsoePXoEERCpXLmyLF++XBo0aBB02XgusHXrVvm///s/KS4uNiyGNn/Dhg2G85mBAAL2FaDNsm/dULLoCRCURc+WLSOAQGoLEJSldv1z9AggEIFAsgZly5YtkwsvvDBoT7KsrCyZN2+etG7dOgLF2K2qHdOiRYsMd3jttdfKlClTDOczAwEE7ClAm2XPeqFU0RcgKIu+MXtAAIHUFCAoS81656gRQMACgWQMyrRbD5s0aSI7duwIKORwOOSNN96w3cD9gQo9aNAgGTJkiOEiEydOlJtvvtlwPjMQQMB+ArRZtFn2OytjVyKCssit9+xxyZa/S2XXLpcceqhDjjjcIYdVS5O0tMi2vXFTqWj/aa+zG2VIenpk22Pt8AS2bC2VHTtdcurJ6aL+dOWFgGkBgjLTVCyIAAIIeAskW1C2a9cuadasmfz666/eB6rzqW/fvjJ69GidOfadNGfOnLKndxqVUHuq54knnmg0m+kIIGAzAdos2iybnZIxLw5BmXnyv7eVyvxFxTJnSYn8/KdLft6l1t2n/ivR2YYWqGSLHFFR5NgqInWOckirZmnS6rwMqVsneOKVn++Saq2d4vx3/7a3fpwj1Y+MMHnTKSaTjAU2bymVx58ukAkfukRUXrlpSY7UqkkdGIsxx1eAoMxXhM8IIICASYFkCsq0sbvat28vCxcuDHr0p5xyiqxYsUJycnKCLmunBbZv3y5HHnmkbpFq1aolmzZt0p3HRAQQsJ8AbRZtlv3OytiXiKAssPmnnxfL7PnF8s6yUvljS+BlTc/NFenU0CFtVHDWunlGWU8l33XvG5QvI6ft702mzVs7L1uOP9Y/YDv83DzZoQV2UXr16eiQMU+qAkf4WvFdsZzVrTDCrexf/fLGDpkxKfIyGRVmq+od+MSoAhn3gQrIPEJQgjIjMaYbCRCUGckwHQEEEAgikExBWe/evWXChAlBjlhEG5fsyy+/lIYNGwZd1o4LaAP2b9y40a9oXbt2lbfeestvOhMQQMCeArRZtFn2PDNjWyqCMn3v2R8Wyv3PFssvsXg+j+p1dlsLh5zdMK3s1r4Z80rlwxUqpPF4rZyZLfVP8Q/KHI3zRPZ4LGjx29vaO+T5UZGHUl+vKJbG11oTlJ12rMgP8ypYfKQi27aXytBnCmX0eyqg9AjI3DsiKHNL8NOsAEGZWSmWQwABBHwEkiUoC/SHts8hy5NPPikDBgzwnZwwn8844wz5/vvv/cqrhYQ9e/b0m84EBBCwnwBtlpR9sUGbZb9zM9YlCvRvYdq0adKlS5dYFylu+ytV+cg77xZK/7HFsmFr3Iqhu+Mvp2RJ47My/OZNmVooX/9QKj+vLZWPflKzyzuh+S0bdIK6q7BpHZFTj3PIUTUcUq2KQ5qdna6736Db8llg585SmTKtSDZscsl69d+MH1QQqN22avalynaBejh6/Tpp0qJpmnTqmGV2zaDLbf+nVIaPLpSnZ+oHZO4NEJS5JfhpVoCgzKwUyyGAAAI+AskQlGm3WrZr105KSnS+fvM53hYtWsiSJUvUALeJO8ZDmzZtyo7B59Bk1apVUq9ePd/JfEYAAZsJ0GbtrxDaLJudmHEqDkHZfvgfVhZLu96FsmV7CBWhRo+4pbVDzjo9TY6ulSZH1VQBk/p5aEWH/L3dJVvVIPAbN6sAa3mJjP9YBUO7Q9i2z6JLXsyUVi0yfaZ6f/x1TYl061cgK9Z6Tzfz6UJ1K+iUsdly5BFwc+dBAAAiLklEQVSx+fvMpTg+/qRIWt1bFNSlaV2RV5/OlpNO9O9RZ+bYjJbRwrvhYwplxHQVkBUbLVU+naCs3IJ35gQIysw5sRQCCCDgJ5DoQdnevXtFG29M71ZE34PNzc0tC5OOOeYY31kJ9blz584yc+ZMrzJXr15d/UFss6+fvUrIBwQQ0ARos/afB7RZ/HtwCxCUiYx9oUD6jlZf9pnpjaWyqhtaOaRb5wxp1TwzpCdRrlpdIq+8USQjZwTuueSuG8+fH4zOlEsuDhyUactrAVSPvk55daF6E8Jr7ze5UlEFfLF+/bmhVI5tm29of8HpIgumVrD0aZPa00lHjCmQYdoYcCYCMrcJQZlbgp9mBQjKzEqxHAIIIOAjkOhB2e233y7jx4/3OSr9j9rtltptl4n+uummm+SVV17xOgzt1hTtFhVeCCBgbwHarP31Q5tl7/M0lqVL5aBs578u6XZ7vsz91lyodHVLh4x+IluOODyyXlfa0xQfGV4gL81X+zW3a3l7aIZcebm52w3Vs5XktHZ5suovk2eSOhzXL9aP+WVy73JcqzxZr/ddozrc7epJk4cfFpm3uxy7d7tk5LgCeXyqCshUR7ZQXwRloYqxPEEZ5wACCCAQpkAiB2WffvqpNG/eXH17GfyvvCpVqsi6deukatWqYUrZZzW9C+1x48aJNp0XAgjYV4A2q7xuaLPKLVL9XaoGZT/+VCynX68GlzcxTlau+tPlg2GZ0qZV8B5doZxP3/1QLGfepMrwX/C1Jj2ULjd0yw6+4IElfvqlRBpcUWDYU8t3Q9+8lSWNzvAfA813uWh8rtU8T/eW11F90qXf7eaP2ahse/fuD8gee9MnINPyN+1uTpOhGUGZkTDTjQQIyoxkmI4AAggEEUjUoKywsLDsqZXaGDdmXkOHDpX777/fzKK2X+aWW26Rl156yaucK1eulPr163tN4wMCCNhHgDaLNss+Z6O9SpKKQdm69SVyQicVIjmD18X5p4nMmZwr2dnRuS3x722lcvH1+fLDH4HLMu6edLn9ltBCoz4DnDLuveBfZmp7Hnh1mgx7RA26FuOX9gCF9PrqyZ3qp9dLPWjT+Xmu5ORE7q498KD74AP3WKonjN7exiGXt8+QFudlSobKBufML5IOd6u0LMhQuwRlXjXEBxMCBGUmkFgEAQQQ0BNI1KBM643Qp08fvUPym1ajRg1Zu3atVKgQv279foWKYELXrl3l7bffPriFww8/XLZt26bGz4j8j7mDG+UNAghYKkCbRZtl6QmVRBtLtaBMC6aO6pgvJSYG1m+jnrI4941cycqK7u/3wkKXXHlLvrz3pXGoNaJXmtx3Z2hBljZY/WEt1fhfZnpMqQCp5KsK6mFLsT25tYcoNOyietX5vJ64JU0euCe04/XZxMGPcxcWyZhXiqX3delqnLcs3WO8rrdTJi829tc2RlB2kJQ3JgUIykxCsRgCCCDgK5CIQVleXp6ccMIJpgevf+6556R3796+h56wn88//3xZvHjxwfJffvnlMmPGjIOfeYMAAvYSoM2izbLXGWmv0qRSULZnj0tOv9SpPx6WT7W0Vp3EtZAsWj3JfHYnWljWvLNTvvrNd87+z4/ckCaDB4YeHA1+Ml8efdW3u5b+Pt5/NlM6tLX29lL9PZVPfeCxfBmm3RLp+VK3Q+5YniPVqsUutXvrnUK55uEDvc48y+LxnqDMA4O3pgQIykwxsRACCCDgL5CIQdmIESNEG5jfzKtmzZry559/SmZmbP/wMlO2cJepW7eu/PZb+V+yo0ePlr59+4a7OdZDAIEoC9Bm0WZF+RRL6M2nSlCm3eJ3zv/y5Ms1wavr6Ooivy+IXUjmLpE22PzRbZ3y37/uKeU/770yTZ56LPSgTAsHK5+r7jE10atMe8Lkwrdj2/u/8jl5smdn+XFq766/wCGvjlP3XsbwtXhpkZzfMzASQVkMKyRJdkVQliQVyWEggEDsBRItKNuzZ48cf/zxsmPHDlNYDz74oDz++OOmlk2UhY466ijZvHnzweKuXr1atPCMFwII2E+ANkuENst+56WdSpQqQdnzrxRIrxFBBqHSKkbdZfnVlCw5u1F8BrZf/mmRtLjJP7Dp1cEh40eGFx4NGpovQ1736bVlcBL+9mG21DleG+E++q+yBypc4X/bZTweLPDVN8XSpJt/WTwVCMo8NXhvRoCgzIwSyyCAAAI6AokWlA0fPlwGDhyocyT+k9LUQBfa2GTHHnus/8wEnqJdVIwaNUqcTqdccskl8vLLLyfw0VB0BJJbgDZLhDYruc/xSI8uFYKynf+65LBWqleVGr8/2OuOyxwydnh4gVSwbZudf/OdTnl5vvd4Wded75DXnguvXGVjlTVXY5WZyAl7d3TIc0+Ftx+zx+de7q4H8mX0TO8A77iaIuuWxLZXm1YegjJ3rfDTSgGCMis12RYCCKSUQCIFZS6XS+rUqSPr1q0zVUdt27aVuXPnmlqWhRBAAAGrBWizrBZle8kokApB2W13O+XFD72DJ926PFTkv6W5csgh0R28X3ffHhO3/l0qNVurYMsjQ+p8jkOmvxJ+gHVrP6dMnGvCQHWk2/VprlSuHF2DggKX5DT2Dy/DebqnB13YbwnKwqZjxQACBGUBcJiFAAIIBBJIpKBswYIFcvHFFwc6HK95M2fOlE6dOnlN4wMCCCAQKwHarFhJs59EFkj2oOz7H4vljKvULXUmMqJXHkqXG7tl26I6r+3plDeXlhe67RkOmftW+EHZmt9L5KQOJrrUqaMfdluaDOwX+nhoocC9Oa1Qrh3kM3i+uuNz9+e5UqlSdEM6vXISlOmpMC1SAYKySAVZHwEEUlYgkYKyzp07ixZ+mXlpg/hv2LBBMjLiM8aHmTKyDAIIJLcAbVZy1y9HZ41AsgdlzTvlySerTFipZw7lfZUrubmxD2n0SlcW8F1ZPmbWeSeLLH83slsS213tlHnflYdvevstm3aIGvv/ywrqbzjDJSKe0bhDnnz9u/dmurV2yOQJ4YeB3lsL7RNBWWheLG1OgKDMnBNLIYAAAn4CiRKUbdmyRWrXri3FxT7f/vkd0f4Jffr0kTFjxhjMZTICCCAQXQHarOj6svXkEUjmoOyvjaVS+wJ1C6OJ123tHfL8qPiENEbFO7pFnmzatn/uaceK/DAvsqBs2SdF0upm/wcF6O1/8uAM6dY1S29WxNN+XVMiJ1/q37vts9eypFmTKKZzAUpOUBYAh1lhCxCUhU3HigggkOoCiRKUPf/889KrVy/T1TV//ny56KKLTC/PgggggICVArRZVmqyrWQWSOagbPgz+TLwBY+BvgJU5A/Ts+S0+vEJaYyK5fm0ymNriPyxNLKgTNvPSefnyZpNRnssn16nlshviyPfX/kWy9/deX++jHnXu16qHyGydXl09le+Z+N3BGXGNswJX4CgLHw71kQAgRQXSJSg7NJLL5X333/fVG1VrFhRduzYIVlZ0fkm0lQhWAgBBFJagDYrpaufgw9BIJmDsmNa5cmGrcEx6h4lsvqj+IU0RiX8eVWJ1O+0v+dVxaoiez+PvIy6Y4MZFGDJi5nSqoW6J9XCl9Eg/qP6pEu/2+M3PhxBmYWVzKYOChCUHaTgDQIIIBCaQCIEZQUFBXLYYYfJvn37TB3c5ZdfLjNmzDC1LAshgAACVgvQZlktyvaSWSBZgzLPkClY/d3TJU1GDonu4PXBymA0X3ta5c7dIqef7JCH74u8jNoIGplN89TjPY32WD79woYOWTDV2ttRdYO6NJEd6kmb1arGb3w4grLyeueddQIEZdZZsiUEEEgxgUQIykJ9ctwrr7wiN954Y4rVJIeLAAJ2EaDNsktNUI5EEEjWoGzA4HwZMdX79j6j+pg+PEM6X5Y6veCfGp0v/SeYs1n9QbbUraMeR2nRS28Q/y7nOWTaS9YGcqEWl6AsVDGWNyNAUGZGiWUQQAABHYFECMruueceGTVqlE7p9Sdt3rxZtKde8kIAAQTiIUCbFQ919pmoAskalDVslyc//GGuVjYuzpGjaqluTSny2rPHJZXPcYqYeD7TLe0c8uIz1oRYRoP4L52YKS2bW3uLZ6hVSVAWqhjLmxEgKDOjxDIIIICAjkAiBGXnn3++LF68WKf0/pOOOeYYWb9+vf8MpiCAAAIxEqDNihE0u0kKgWQNyhynqdsLC01U0SEirhWRj/1lYk+2WqTvQKeMneUKXibVmWzH8hypVi3yIFFvEH+rxl4LfiCBlyAoC+zD3PAECMrCc2MtBBBAQBIhKDvyyCNl+/btpmrrqquukqlTp5paloUQQACBaAjQZkVDlW0mq0AyBmVb/y6Vmi3zTVXZJWc55IMp1vSYMrVDmyy0fkOJHHexelCAiazs0R5pMqh/ZOOjGQ3iP7xnmvS/K7JtW0FKUGaFItvwFSAo8xXhMwIIIGBSwO5B2d9//y01aqhnkpt8abdo9uvXz+TSLGaFQFFRkTz44IOyZs2asocuPPTQQ3LcccdZsemobSMRyxw1DDZsqQBtlqWcUdmYnf/9az2ily1bJhs3bhTtXNq2bZs4nU6pVKmSVKlSRapWrSonn3yynHnmmXLiiSdKWlrkvWyighzCRpMxKFu8tEjO71lkSqFvpzQZPSz+QY2pwlq8UKcbnDLrCxNJmepwV/BFrnqaefiD7b/xdqF0e8TnXk/1z+fvZTly5BHx/3dEUGbxycXmygQIyjgREEAAgTAF7B6ULVq0SC688ELTR/fZZ59Js2bNTC/PgpEJ7NixQ6644gpZunTpwQ09/PDD8thjjx38bLc3iVhmuxlSHmMB2ixjGzvMseO//99//13GjRsnc+fOLfvCwaxT5cqVpVu3bnLHHXdIvXr1zK5mu+WSMSgb83yB3PlsiSnrwTemySMDUjMo+/LrYmna3cz9qSKTHkqXG7plmzLVW0hvEP/Lmjhk1mv26M1HUKZXa0yLVICgLFJB1kcAgZQVsHtQNn78eLn99ttN1Y/D4ZB9+/ZJbq49/ugxVegEXcjlcon2dNEHHnigrMeD52Hceeed8uyzz3pOssX7RCyzLeAoREgCtFkhccVsYTv++9d64Q4ZMkTeeustKSkxF6rogWm/+7QvlPr06SPt27dPuF5myRiU3drPKRPnmugppSp0TL906XNb+AGQ3jmRSNPMPvSgtrq54M+l4Y3lZjSI//zxmXJRm/gO4u+uK4IytwQ/rRQgKLNSk20hgEBKCdg9KHviiSdEu5XPzIuB/M0oRbbMhg0bysaAmzx5svz000+6G9Mu1saMGaM7Lx4TE7HM8XBin9YI0GZZ42jVVuz673/ixImifamg3Vap96pYsaLUrl1bsrOzy27B1G7DNBOmHX/88TJhwgS56KKL9DZry2nJGJS1v8Ypc781F5RNeTRDrr0qy5Z1E4tCzZxdKJ37+9wSabDjhRMy5YLWoQdbeoP4SyWR0i8riMqZbfEiKLNFNSRdIQjKkq5KOSAEEIiVgN2Dsv79+8tTTz1likP7Rn3BggWmlo3mQv/884/89ddfovVgSE9Pl4yMDMnMzCy74MnKyir7qV38aO+16WZfpaWlZRdV2oWV9l9+fr4UFhaKNt6O9qpbt65UqBDet61GZdAeoqDdFrRixYqygEy7tVU7rkCveAdliVjmQJ7MSywB2qzy+qLNKrdwv9Pa6xtvvFHeeOMN96SDPxs3biw33HCDdO7cWbQHQni+tDZ//vz5MmPGDHnnnXekoEANgm7w0sYt0wLbgQMHGixhr8nJGJS1uDxPlv9iznnOmExpf5H5vwXMbTVxllJ/2khWszwp2R28zK3riyyeHtrfOUaD+D92U5o8fJ99bnklKAte/ywRugBBWehmrIEAAgiUCdg9KLvtttvkxRdfNFVb2i2a2jgv8XxpFzFXXnmlaBeIZl7aLTNaaOYZnGnhmTa9uLi4LATTAjHtIkkLxQK96tSpUza+jbZuOC+t5963335bFsDt2bOnLCDbvdvEX64+O4tlUJaIZfbh4mOSCdBmFZW1IbRZ/ie29iVD9+7d/UKypk2bitbDrH59lQKYeK1du7asN9qcOXMCLt2lSxd58803y76sCbhgnGcmY1B2Zrs8+e4Pc7CLX8iU1i1TNyjTlEIZ0+3nWdlySr10c7hqKd1B/NWfSZuX5EjNGvEfxN99IARlbgl+WiqgfvHwQgABBBAIQ0D9oa11D9L9b9iwYWFs0dpVrrrqKt2y6ZVZPfHS2p2HsbXLLrvMdHn1jiHSaf/++28Ypd6/SsuWLU2XXYVxhsuqoCzsMoS6YiKWOdRjZPnEEqDN0v99YtS2pVKbde+99/q1m+r2S5f6EiSsk3zo0KF+2/N1vu+++8LadixXUrfqGx7HtGnTYlkUy/Z10vn7XHKSuf/efb/Asv0m6ob++6/UJQ3MeV1/e15Ih3n2Jf7bbds1tG2EtMMwF/7y66Kg58ymzSVhbp3VUlVAuw2EFwIIIIBAGAJ2D8o6duxo+Ae07wWB+uY8DAFrV2nRooXp8vqWP9LPVatWjehggoVOtWrVcvXq1cv19ddfu9RA/obHaaegzI5ljqiSWNn2ArRZ5oOyVGqz1G3rLs8vGLT3U6dOjfh8HjlypGFb7P6d8u6770a8n2huIBmDsiYd/cMZo+Dspdfyo8mbMNvu/4gzaFBUZnjKPtff28wFRqtWF+tu84N54YXT0cQkKIumbupum1sv1W9CXggggEA4Ana/9fKaa64peyKYmWNbvHixtG7d2syiUVvmvffekxEjRog2TpY2Fo02+LJ2C6V2G5L2X6BxZUItlHaLZqVKlcpu26xRo4aongPStWvXUDdzcPl+/frJypUry8ZNq1KlStkYOdp2GzRoII0aNZKaNWseXHbSpEnSo0ePg58938Ty1stELLOnFe+TT4A2y7hOU7XN0m7FP+uss+S77747iKO+dBDtCamRvtTln2hjm33zzTeGm6pcubL89ttvcsQRRxguE88ZyXjrZacbnDLrCy3DDP4a3jNN+t9ln7Gygpc4Okts2lwqR1+Qr0bYD779h69Lk8ceCG6mO4h/RZGSryqop8MG308sl+DWy1hqp9C+Ujcj5MgRQACByATs3qNM68Gkfp2Z+k+FPJFhxGBtFZa5Nm7c6Fq6dKlL6wlwwgknmDo2t0Hz5s1dc+fOdakHBsSgtMa7sEuPMuMS+s9JxDL7HwVT7C5Am+XdXtNmuVyzZs3yaufVE5pde/futexU/vjjj7227/594fnTDkMpGB1wMvYou61fnm5PJr1eZfcNchrRpNz0q2816Xb6Pld+fmlAH22+nObfs++hx+3pTY+ygNXJzDAFuPUyTDhWQwABBOwelN1///1BLwDcFwPr1q1LuAr96aefTB9f7dq1Xbt27bLFMSZi6JSIZbZFZVOIkARos8qDMtqs/aeO77h1quezS/Xmdu3cuTOkcyvQwldccUXA3yVaOKd6OAfaRNzmJWNQpoUxeqGY3rQb+9hvvKx4nQzf/RB8nC634QuTAt+yOmVqgW4d/PmXPf8dEJTF66xL7v0SlCV3/XJ0CCAQRQG7B2XDhw8P+Me/OyTTfq5atSqKUtHZtLodxvTxvf3229EpRBhbTcTQKRHLHEbVsEqcBWizyoMy2iyX6vWS7zrkkEN02/kjjzzStWTJEkvOWDO/Sz744ANL9mX1RpIxKBvzfL5uSOMOeTx/duhGUOZ5Tp1zqX8vME8v9/vq5+3zXM3vvd4g/m2uCLyO30ZiOIGgLIbYKbQrgrIUqmwOFQEErBWwe1CmxsLSvcDwDMjc75ctW2YtTgy2NnPmTFPHl5OT49q3zz5/4CVi6JSIZY7BKcguLBagzdoflNFm7T+xVq9eHbCNP/TQQ12RPPnT8/TVgjf370O9n4MGDfJc3DbvkzEoe+sd/d5M7pDH8+cJre3zu90OJ8Wc+YWmQ8YPF+gPym80iP/M2fZ9wihBmR3OvuQrA0FZ8tUpR4QAAjESsHtQpgY/DviHv+fFwDvvvBMjNet2079/f1PHp92qY6dXIoZOiVhmO9U5ZTEnQJu1Pyijzdp/vixatChoG//888+bO7mCLNWqVauA++revXuQLcRndjIGZT/9ov+0Rc+A7OD7evtce/YEHm8rPjUTn72WKorDzjHXq+y8/+mHjH0H6tz62nCfq6goPsdkZq8EZWaUWCZUAYKyUMVYHgEEEDggYPegrLCw0JWdnR3wj393WDZu3LiEq1dtoGt3+QP9fOyxx2x1bIkYOiVimW1V6RTGlABt1v6gjDZr/+nyxhtvBG3jBwwYYOrcCrZQ7969A+7rvPPOC7aJuMxPxqBMgxQ14PzBMOykwO/nLtTvGRWXCrHBTrXxx8za/fiTd/pVNoh/A3/v/o/YcxB/NzdBmVuCn1YKEJRZqcm2EEAgpQTsHpRpldGoUaOAf/y7A6aHH344oepOG7vGbAi4fPlyWx1bIoZOiVhmW1U6hTEtQJslLtqs/aeLNmi/+3eU0c/HH3/c9LkVaMFAgZO2b+3hCnZ8BSr3tGnT7FhkU2W67DqTT3BUIdqAwfYOcUwdsIULOZ3qiZUmg8Zre3qP8WY0iP/aP4otLKH1myIos96ULarAHgQEEEAAgfAEEiEou/XWW4NeaGgXAbfcckt4CHFaSxtTzejCyXN6hQoVXFovFTu9EjF0SsQy26nOKYt5Ados2iz32fLnn3+6HA5HwLZ+3rx57sUj+vnyyy8H3E+dOnUi2n60Vk7WoGzsC+Z7RZ3ZTv8WwmiZJ8J2H35C5/ZJvZ55J+9zbdla/iTLRu39e5Ode5n9fQnKEuGsTLwyEpQlXp1RYgQQsIlAIgRls2fPDvjHvztUOuWUU2yiaq4Yjz76qKnjuuiii8xtMIZLJWLolIhljmGVsisLBWizaLM8T6du3boZtvVNmjSx7IuQJ5980nA/2u/Jdu3aeRbLNu+TNSj7eVUI45SpsEe7ZZBXucDf20pccop/6KV3S+b9j+7vkffLr/rmU6fbdxB/9xETlLkl+GmlAEGZlZpsCwEEUkogEYKygoICV5UqVQJeALjDsnXr1iVM/bVs2dLUMQ0bNsx2x5SIoVMiltl2FU+BTAnQZtFmeZ4oxcXFrn79+pX1LNOeBur+fdW+fXvXxo0bPReN6P2NN954cNvufXj+vPPOOyPafrRWTtagTPOSM8wFPVr4oz3tkZe3wA13mLx99bR9rry8UlefATrLq3k265TvfZAHPhGU6bIwMUIBh7a++kXACwEEEEAgRAE1fopccsklumupgEYGDhyoOy/WE2+66SZRQUfQ3ao/uKVPnz5Bl4v3Anl5eVK1alVRt1QGLcoXX3whqtdB0OViucCkSZOkR48eurvU/LV6sNsrEctsN0PKY16ANos2y/ds2bx5s6xdu1bS0tKkevXqom6F9F0kos/HH3+8/PHHH4bbeO6550QN+G84P14zxo4dK3379tXdvRqjTLp06aI7LxEmDhqaL0NeLzVV1HNOEvn0vQqmlk2VhVatLpFTLiswdbhP35Eu97xYIuLzZ9VdndPkmSdyTG0jngt99U2xNOnmU3ifAm1akiO1aqb5TOUjAsYCBGXGNsxBAAEEAgokSlC2cOFCUbcgBjwWbaa2zPz584MuF+8FtDK2bds2aDEqVqwo//77r2RkZARdNpYLJGLolIhljmWdsi9rBWizaLOsPaMCb2316tVSr169gAt9/vnn0rRp04DLxGNmMgdl2/8plSNb5ouo/MbM6+s3s+SsM+31b8dMuaO5zAVX5slHP4a/h1/fz5aTTkwPfwMxWpOgLEbQKbYbgrIUq3AOFwEErBNIlKBM6zh8+umny8qVKwMefFZWlmzfvl0qVaoUcLl4z7zuuutk8uTJQYuhhWlz584NulysF0jE0CkRyxzremV/1gnQZllnacWWkv3ff8+ePeWFF14wpDr77LPlq6++MpwfzxnJHJRprj36OmXSAnM3P3U42yHvT86NZ3XYbt+LlxbJ+T2LwipXoxNEvpmTGL30CMrCqmJWCiag3XrJCwEEEEAgdIFEGKPMfVTvvvtuwPFX1O+KsvkPPfSQexVb/tyxY4fLc5wad7n1ftpxfDINNRHH+0rEMtvyBKZQpgVos0xTRX3BZP73v3Xr1qC/U9QtjFE3DncHyTxGmWayarX+APN6g9Jr01b/VhwuZdKud2xL82O9ebq+/qb9B/F3VxpjlLkl+GmlAD3KgiWJzEcAAQQMBBKlR5m7+GeddZasWLHC/VH3Z25urqxZs0aOPvpo3fnxnvjMM8/I3XffbaoYn332mTRr1szUsrFcKBF7ZyRimWNZp+wrOgK0WdFxDXWryfzv/8EHH5ShQ4cakmhjl2m/E9PT7Xn7WbL3KNMqpt3VTpn3nbleZVec65B3XrZXr7K1f5TIzp0uObtRfG4LnfxWoVz3aLHhOa47I0vE+VWu5OQ4dGfbbSI9yuxWI8lRHoKy5KhHw6PQBj/VbrfSfu7evVt27dpV9nPPnj2SmZlZNiB2tWrV5PDDDy+7NatBgwZl0w03GMaMvXv3ypAhQ7QnrMqgQYPk0EMPDWMrrIKA/QQSLSgzO7ZX9+7d5fXXX7cduNaGnHzyyaKNJxPsVaFChbL2Tmvn7PZKxIvORCyz3eqd8oQuQJsVulk01kjWf//aGJZaEKb9baz3cjjUrXzvv2/40B69dWI9LRWCsr82lkrtdmqsMpN3EI7ply59bsuOdVXo7q/s1sc+quBqnLXZT2dIx3YqgYrxq0jtPqtJnoj6n9lX744Oee4pewWOgcpOUBZIh3lhC1jZPY1txVdABVIu7VYF7RHWrVu3dh122GGmbrVSJ8/B5bRbms455xzXyJEjXSpci/iAVDjnUk+cO7j98ePHR7xNNoCAXQQS6dZLt5l6mtzBf4+e//Y936uLA5caj8W9im1+jhgxImjZ3cdxwQUX2KbcvgVJxNuYErHMvu58TkwB2qz411uy/vvv2rVrwN8pAwYMiD9+kBIk+62X7sN/e0aBy/O2wIDv6+1zLVxc6F41bj8nvJzvElWWg2U9eZ+rME7FGvq0s7wc6hbVg2UyeL/y58S6hZVbL+N2mif1jrVePrwSWEA9yto1evRol3panSs7OzvgL3z3BaTZn+oR3K4OHTq41K1aYQnl5+e7WrVq5VWmqVOnhrUtVkLAjgKJGJT9999/rhNPPNHr36Vem1C7dm3XX3/9ZRv2pUuXutStL0HL7T4W1YvVNmX3LUgiXnQmYpl93fmcmAK0WfGvt2T89z9lypSAv09atGjhKi62f1iQKkGZ9q+gW8+8oAHPwQCofvzGKyspcbn6DPAv6+gJ+XH7x/zvv6UuOTV4QKb5nXrhvriVM9wdE5SFK8d6gQQIygLp2HReUVGRa+bMmWXhmNbzw31hGM2fXbp0camn4ZkWKVG/JTp37uxXNit6qZkuBAsiEGWBRAzKNBKtt5i6JdHv36dvG3LKKae4tMHz4/3atGmTq3r16kHL61n+jz/+ON7FNtx/Il50JmKZDSuAGQknQJsV3ypLtn//f/75p6ty5cqGv1Pq1KljyV0Vsai1VArK9u4tdWU1MRf2aIFPxab7XOvWxzbs3LGz1HV+F/8yjnjWGYvTIeA+et3jH94dDBY9epa9/Hr8Ar2ABxBg5tLlhUFD1E2bVYLJC4EQBP4fAAD//4RXk9IAAEAASURBVOydB5hTRdfHzy67yy7S+4LSm4hKV3rvvCh9UUBp0l5QXkAQAVEQFBAE6U0FBAFpAtK7oCAICop0pPeysL3cbyb7ZUm5yb1J7k1u+c+z++SWKWd+J5kk/5yZCRJYIiRdELh79y7NmDGD5s+fT9evX/fI5qCgICpTpgwVKVKEcuXKZfnPmjUrPX78mG7cuGH5v3DhAl29etVlvc8++yytXLmSqlWr5jKP9Ubfvn1pzpw51lPLY6lSpej06dN213ACAnom8NNPP1GLFi1EuzBhwgQaPny46D0tXFy6dCm99dZblJqa6tYc/nrfsWMHZcqUyW0+tW6eP3+eXn/9dTp58qTsJsLDw+nRo0cUFhYmu4w/M3799dfUvXt30SYHDBhA06dPF70XyIt6tDmQvNC28gQwZinPVG6NRnr9x8XFUd26denw4cOi3S9dujTt2rWLChQoIHpfaxe/+uorGjhwoKhZ/DN7+/btRe/p9eK5CylUslMC0SOZPQglWjY6hDq1V/fzAP8oNffrBOo3PYWImWebPukRTKOGhtteCsjx+YspVKI5M87dN/8QopjDEewzX1BAbPS20fWbEun1wclui1/dFU4FCwS7zYObIGBHgAtlSNomwL7wCaNHjxayZMnChzZZ/8HBwUK9evWEsWPHCtu2bRMePnwoq5PHjh0TRo4cKZQtW1a0ndDQUOHLL790W9dHH30kWrZXr15uy+EmCOiNwKZNm0Sf6/x1yoQyzXdn3rx5AhPRXfbBOt5UrFhROHXqlN/7s3HjRiF79uxO9pUvX96t3Xzs03JatGiRU5+srJlQpknT9WizJkHCKJ8IYMzyCZ/XhY3y+mc/DAlt27Z1Of7yz743b970mlMgCrIfVlz2hwllgTBJ9TZv3koRiteLEai0/P/OfWKFBw9SVbFt38+JQpE64raM+CROlTa9rbTFm7FuufV8N9bbqgNabvL0eLf94s+Vo8eSAmojGtcfAdKfyeaxODY2Vpg4caLAIsBcvglav1xZH6tUqSJMnTpVYBFnPoPiIkDx4sVF2+7Xr59o/bNmzRLNz+1jvwaLlsFFENArAb0LZZz7tGnTXL5mreMKf4yIiBBYRKtfXBUdHS2waDxRMaxhw4YCiwgQevbs6dLuMWPG+MVObxvR45dOPdrsrX9QTtsEMGb53z9Gef0PHTrU5fvGq6++Kty+fdv/cH1s0YxCGUcWHZ0qVG8lLk65FNBeiBEGDo8Trl5L8ZG6IDDNVTjwS5LQ5m3XwtMHH2tLJOOd/vlgkltBSa9iUtMo136wPh9+WJfgs99RgbkIQCjTqL/ZVCehWLFiLt/Qbb/A8ugxFlot8GgwpRP/Qsq/dGbMmNHJlk8//dSuuVWrVgncFlvbbI+vXLlilx8nIKB3AkYQyrgPVqxYIWTOnNnla9f2ddy4cWPh0KFDqrjuyZMnlkg8Vz8OdOrUSUhISPugc+7cOYFHuNraZj3ev3+/KvYpVakev3Tq0Wal/IV6tEcAY5Z/fWKE1/+IESNE3y/4+0bnzp2F+Ph4/0JVqDWzCmUcX0JCqjB8TJxAz3somLH87XvECkuWJwjXrssXzVJY1t17E4U+/2OiTAXXbUbWjLEIUgq5WPFqXmgkbnsJFqWnx3T6bLJAZcT7ZBXJ+GPfwfqMltOjT4xiM4QyjXny/v37Alu7xuWbufWLIH/kolTXrl39MiVq69atlogS2/b58eLFiy0E2XoOomKaNT+PTEMCAaMRMIpQxv3Cp1a+8MILssYe/rpma5dZBLakJN9C2flUmAMHDgjvvfeekCdPHtH22Vpjwrhx49gvuPbTJrhYbx1jrI9c8EtMTNT0U02PXzr1aLOmnwQwzmcCGLN8Rii7Ar2//ocNG+b0XsHfM/jSA+PHj5fNQYsZzSyUWf3x16lkoUJTaaHEVjSxPc5RPUZ4k03NHDIqTpj4ZZzw9ZJ44Zul8cLnU+MsEWhtu8UKVVuy+l+WaKNsjDBhSpzg48cia7dUe1yxOkE0qmzWAv2JxZevpAgFa0n4xWaK7ur1iCpT7YllwIohlGnIqXw9nvz584u+mVu/BFofebTZvn37/Gr9zp07Bbagt519fDrWsmXLBLYxgN11q53WRy7+IYGA0QgYSSjjvuERXb179xYyZMjg9vVsfV3zRz5m8UgvPu2aLbjvJGY5+jw5OVn4448/hAULFgjvvPOOULBgQbdtVa5cWThx4oRjNZZzXtcrr7xiV75Vq1aiebV0UY9fOvVos5Z8DlvUIYAxSx2ujrXq+fU/ZMgQu/cI6/tXtmzZhB9//NGxq7o7h1CW5jL+O9rcr9k6VeXliya2Ypmvxw07xAiX/pUfnRbIJxr76CRQFQdOTOTj01m1nLiPHz9OFU6dThbWbkgQ+HpqxOz2yHcs8uy1rrHCwsXxwu/Hk4QbN1MEzgMJBMQIBPGL7E0DKcAEJk+eTOwXL8kd6LiZ7MslffHFF8QiJ/xuNdsYgJo2bcoFVo/a/vbbb4lFv3lUBplBQOsE9LzrpTu2fIdJtpYLbdmyxV020Xt8N93IyEjLzrq5c+cmtgkJsemSdOfOHWLTry0767KIL9GythdffPFFy5gYFRVFTLizvWV3zHfFrFGjBt26dctyfe7cuZYx0i6Txk6+xq6XGvMIzNE7AYxZ6npQj2MWJzJ48GCaMmWKExwWPU1r166lkiVLOt3T2wWz7Xop5Z+kJKLV6xNp0qJk+v2CVG4f77MdNTu8GkTdo0KoSQN2oqM0ZUYCDZ7Bduj8/9SlfhAtnhVhPdXcY6OOsbTjDxXNYpt89moaRPOmapeBir1H1a4IiKlnuOY/AnyKUI8ePUR/7WI+s7vOtqoWNm/e7D/jXLTUt29fO7sc7RQ7v3TpkovacBkE9EvAaBFljp7gO+Y2a9bM5VpgYq91X67x6eR16tQRNmzYIBmZZmsrjzhjYpnAd+fk09e1nvQYnaFHm7X+PIB9yhPAmKU8U16jHl//7777ruhnVb6m7+PHj9UBFYBaEVHmGvqfJ5OE3oNihYhXPYw6spmq5xStVC5GiOoVK6zflGBZI81169q+wyOziPXF2r9fD/u2jIbave3QU3qxfmtfvH3s0hdrmKntR73Vj4gy9q0uUOnBgwfEtqmm3bt3S5rApjbRzJkzKUeOHJJ51c7APmAQ/zWOR4fISYULFyYmlMnJijwgoCsCRo0oc3TCw4cPaf369cQ27KDt27eTnIgwxzpcnbPpL9SkSRNq2bIlMVGOeBSa0ZMeozP0aLPRn0fon2sCGLNcs/Hmjp5e/+yLGA0cOJDYLs12XWXrkdHYsWPpww8/tLuu9xNElMnz4M1bqbT352TafSCFdh4X6NxDVi6e/ae6KB9CFJGF6EX2kaRQ3iAqFBlEdWsEWyLHwsJY+JEB0pBR8fTFD6lUKB/Rv3syGaBH6AIIKEsAQpmyPGXXdvbsWcsXwzNnzrgtw3abJP4BhQtlWkrr1q2j1q1byzLprbfeom+++UZWXmQCAT0RMItQZusT/gWUreti+efTHq9fv26ZVsm/nEgltig/sTXNiEXHUvXq1S1jYM2aNYntXilV1FD39fSl0wpejzZbbcejuQlgzPLd/3p5/fP3of79+9Ps2bPtOs2XAFi6dCmxNSztrhvhBEKZb16Mjhbo9t1UuntXYMtECJQ3bzAVjAymrFmNIYb5RgelQcDcBCCUBcD/e/bssUSSsSlCblvnb+z8C2ndunXd5gvETbb7HJUoUYIuXrwo2TwL2adu3bpJ5kMGENAbATMKZWI+Yjtf0o0bNyyiGRfO+D9fl4ztYmlZr4yvWcb/c+bMSfxXfbMnvXzptPWTHm22tR/HIGBLAGOWLQ3pYz28/rlI1qdPH5o3b55dh4oWLUpsOr9lJoTdDRknvN9jxoyx5GzcuDHNnz9fRin/ZoFQ5l/eaA0EQMBEBNgbC5IfCSxcuFDWej/sC6Zw5MgRP1rmeVOffvqp6PoP7OVjd51FnXheOUqAgA4IGH2NMh24QJcm6nG9Hz3arMsnB4wGAQ0S0Prrn/14K/Ts2dPusyf/LFqtWjWBbfTiNdF+/fql18mioL2uR82CWKNMTbqoGwRAwMwE+O6FSH4i4O6Dhq24xLes/vPPP/1klffNsAgSISQkJP1DhG0frMfPPfec9w2gJAhonACEMo07SKPmuXsvGDBggCat1qPNmgQJo0BAhwS0/PpPSUkR2KwFp8+iHTt2FOLi4nyizcUx6+dZLsRpMUEo06JXYBMIgIARCEAo85MX2XRLWZFkbA0fYdeuXX6yyvdmGjZsmP4hwvphwvbxzTff9L0R1AACGiUAoUyjjtG4WVr+0ukKnR5tdtUXXAcBEPCMgFZf/1wk69q1q9Pn0JEjR3q0c7IYDTY9V4iIiEive8qUKWLZAn4NQlnAXQADQAAEDEoAQpkfHHvu3DkhV65c6W+2tkKS7TFbu0dYtmyZHyxSrolRo0a57Rdbz0G5xlATCGiMAIQyjTlEJ+Zo9UunO3x6tNldf3APBEBAPgEtvv6Tk5OFzp07230G5bMcuK1KpIMHD9rVvWXLFiWqVbwOCGWKI0WFIAACIGAhAKFM5SfCgwcPhDJlyti92dqKY7bHbDtrla1RvvqNGze67Rvb1VP5RlEjCGiEAIQyjThCZ2Zo8UunFEI92izVJ9wHARCQR0Brr38ukrHd4O0+f2bNmlXYtm2bvA7JyDVixAi7+i9fviyjlP+zQCjzP3O0CAIgYA4C2PWSKVVqJfZGTs2bN6ft27dLNlG2bFk6evQohYeHS+bVUoY7d+6wrZTzippUoEABunbtmug9XAQBIxDArpdG8KL/+6CHHeQcqejRZsc+4BwEQMA7Alp6/fPP1iySjFasWJHeGf55k78fv/zyy+nXfDl4/PgxlS5d2rKbM6+H70IfHR3tS5WqlcWul6qhRcUgAAImJwChTMUnANsth2bPni3ZAluXjA4dOkTly5eXzKvFDGzBfrp69aqTaVFRUbR8+XKn67gAAkYhAKHMKJ70bz+09KVTbs/1aLPcviEfCICAewJaef1zkYxFktEPP/xgZzDbUZ4KFSpkd83dCYuFsNxmu2USW4uMEhMTKSYmhh49ekR3796lHTt20IEDB9KrqFq1quVzevoFDR1AKNOQM2AKCICAsQiYI3DO/710FwrNnkF24dyfffaZ/w1UsEUm8Nn1x9o/JhIq2AqqAgHtEcDUS+35RA8WzZgxQ3TM5GNnjx49NNkFPdqsSZAwCgR0SEALr38mZglt2rRxOXZaP3uq8fj2229r1mvuvm+sXLlSs3bDMBAAARDQOgFElLF3VKUTn2rZrFkzYrvxSFZdu3Zt2r17NwUHB0vm1WqG+vXrW/rgaN+pU6eIrc/meBnnIGAYAogoM4wr/dqRoUOH0uTJk0XbbNCggSWaQfRmAC/q0eYA4kLTIGAoAoF+/fOIrw4dOtD69esDwnXSpEk0ZMiQgLQt1SgiyqQI4T4IgAAIeEcAQpl33FyW4usa8PXGxKYiOhZi204TF5MKFy7seEtX523btqU1a9bY2ZwvXz66efOm3TWcgIDRCEAoM5pH1e/Pzp07qUWLFpSQkOCysZkzZ1Lfvn2J7YTsMo8/b+jRZn/yQVsgYGQCgX79c5GsXbt2tGHDhoBhZtHjljWHA2aAm4YhlLmBg1sgAAIg4AMBCGU+wBMr2r9/f5o1a5bYLadrw4YNIzbt0um63i6wqULEdkSyM7t9+/bEQr7truEEBIxGAEKZ0TzqfX9u375tN+bxtXT42jdcELt3755lY5Pjx4/T2bNnZTUSGRlJlStXtvyQkjt3buI/rISGhlJISEi6gMajefkPM94mPdrsbV9RDgRAwJ6AHl7/bFqORSRz/DHWvifqn128eJGKFCmifkNetAChzAtoKAICIAACMghAKJMBSW4WvvBnrVq1iL+xS6Xs2bPThQsXKEeOHFJZNX9fTBxk61kQv44EAkYmAKHMyN71rG9HjhyhKlWqeFbIx9zz58+nnj17el2LHm32urMoCAIgYEdAD69/vqh+o0aN7Oz298kzzzxDfLaIViJ8HfsPocyRCM5BAARAQBkCEMqU4WjZMYfvWsmnUspJ48ePpw8++EBOVs3n6dWrFy1YsMDOzhMnTlC5cuXsruEEBIxGAEKZ0TzqfX/08KXTsXd6tNmxDzgHARDwjoAeXv985/Q33njDuw4qVKpSpUrEWWk1QSjTqmdgFwiAgN4JQChTyIM8gmrAgAGyasufPz+dP3+eMmXKJCu/1jNFRUXRihUr0s3k04R4SL9Wf31LNxQHIOAjAQhlPgI0UPHr169ThQoV6MmTJxQeHm75z5gxI4WFhaX/86mT/J9f41Mo+X+GDBksFPjmL9bpmnzKJl+Xx/rIj/k/n8Zp/Y+Li6O5c+f6FFGmR5sN9JRBV0AgoAT08Prnm2M1btw4oJy6dOlCixcvDqgN7hqHUOaODu6BAAiAgPcEIJR5zy69ZGxsLBUvXlz24vV8oeZ+/fqll9f7Ad+lbdeuXendYNt30+rVq9PPcQACRiUAocyonkW/QAAEQAAEQED7BCCUad9HsBAEQECfBCCUKeC3iRMnEl+YX07iCzT/+++/lqgCOfn1kKdUqVJ2C1RPmzaNBg4cqAfTYSMI+EQAQplP+FAYBEAABEAABEDABwIQynyAh6IgAAIg4IYAhDI3cOTcio6OpmLFill2NZOT/8MPP6Rx48bJyaqbPAULFiQewm9Np0+fJi6eIYGA0QlAKDO6h9E/EAABEAABENAuAQhl2vUNLAMBENA3AQhlPvrv888/p+HDh8uqJTg42LI2mVa3mJbVCZFM/E16ypQpxNfMadGiBS1cuFAkFy6BgPEIQCgznk/RIxAAARAAARDQCwEIZXrxFOwEARDQGwEIZT54TBAEKlGiBF24cEFWLU2bNqXNmzfLyotMIAAC2icAoUz7PoKFIAACIAACIGBUAhDKjOpZ9AsEQCDQBCCU+eCBbdu2UZMmTWTXsGbNGmrdurXs/MgIAiCgbQIQyrTtH1gHAiAAAiAAAkYmAKHMyN5F30AABAJJAEKZD/Tbtm1LXPySk/gi/pcvX6aQkBA52ZEHBEBABwQglOnASTARBEAABEAABAxKAEKZQR2LboEACAScAIQyL11w48YNKlSoECUnJ8uqYcCAATR9+nRZeZEJBEBAHwQglOnDT7ASBEAABEAABIxIAEKZEb2KPoEACGiBAIQyL70wZ84c6tu3r+zSW7dupcaNG8vOj4wgAALaJwChTPs+goUgAAIgAAIgYFQCEMqM6ln0CwRAINAEIJR56YFWrVrRhg0bZJXOnDkz3bt3j8LCwmTlR6bAEoiOjqapU6fS8ePHKSYmhlJTUwNrkIZaDw0Npfz581PRokVp4MCBlD17dg1Z539TIJT5nzlaBAEQAAEQAAEQSCMAoQzPBBAAARBQhwCEMi+4JiQkUK5cuSwiipzibdq0odWrV8vJijwaIMA3aOAbNSC5J1CzZk3av3+/+0wGvwuhzOAORvdAAARAAARAQMMEIJRp2DkwDQRAQNcEIJR54T5Pd7tctGgRdevWzYuWUMTfBFJSUigiIoKSkpL83bTu2gsODiYefffMM8/oznalDIZQphRJ1AMCIAACIAACIOApAQhlnhJDfhAAARCQRwBCmTxOdrkGDx5MU6ZMsbvm7uT69evEd71E0j4BQRCoQIECdPPmTe0bG2ALs2TJQvfv3zf1Tq4QygL8JETzIAACIAACIGBiAhDKTOx8dB0EQEBVAhDKvMDboEED2rVrl6yShQsXpkuXLsnKi0zaIMBF0Pfff594dBmSOIGgoCD66KOPLP/iOcxxFUKZOfyMXoIACIAACICAFglAKNOiV2ATCICAEQhAKPPCi3nz5qU7d+7IKtmxY0f6/vvvZeVFJu0QuHHjBv3777/EH/madEhpBLhAxqdavvTSS1SoUCHTY4FQZvqnAACAAAiAAAiAQMAIQCgLGHo0DAIgYHACEMo8dPCtW7csu/7JLcajkwYNGiQ3O/KBAAjoiACEMh05C6aCAAiAAAiAgMEIQCgzmEPRHRAAAc0QgFDmoSt27NhBjRo1kl3q4MGDVK1aNdn5kREEQEA/BCCU6cdXsBQEQAAEQAAEjEYAQpnRPIr+gAAIaIUAhDIPPTFr1izq37+/rFJ8mlpMTIxlF0VZBZAJBEBAVwQglOnKXTAWBEAABEAABAxFAEKZodyJzoAACGiIAIQyD53x6aef0siRI2WVwkL+sjAhEwjolgCEMt26DoaDAAiAAAiAgO4JQCjTvQvRARAAAY0SgFDmoWP4boiTJk2SVYpP0dy2bZusvGpmunv3Ll25coUEQaAMGTJQSEgIhYaGUsaMGSksLMzyaD3m1+Wm1NRUiouLS/+Pj4+nxMRESkpKslRRqlQpypQpk9zqkA8EdEcAQpnuXAaDQQAEQAAEQMAwBCCUGcaV6AgIgIDGCEAo89AhvXv3pnnz5skqxadozpgxQ1ZetTKtXr2aOnToQFzUkpP4dFEumtkKZ1w849eTk5MtIhgXxLhAxkUxd6lEiRJ05swZS1l3+XAPBPRKAEKZXj0Hu0EABEAABEBA/wQglOnfh+gBCICANglAKPPQL1FRUbRixQpZpbSw4+Xrr79O69evl2WvGpkePHhA2bNnV6Nq1AkCAScAoSzgLoABIAACIAACIGBaAhDKTOt6dBwEQEBlAhDKPATcqlUr2rBhg6xSy5Yto06dOsnKq1amOnXq0L59+9Sq3m29OXLkoPv377vNg5sgoGcCEMr07D3YDgIgAAIgAAL6JgChTN/+g/UgAALaJQChzEPfvPHGG7R8+XJZpXbt2kX16tWTlVetTDyabOLEiXTnzh3LtMmUlBTLFErr2mIJCQmKNc2naGbNmtUybTN//vw0dOhQ4hF4SCBgVAIQyozqWfQLBEAABEAABLRPAEKZ9n0EC0EABPRJAEKZh37r168fzZ49W1apEydOULly5WTlDVQmvt7YvXv36Ny5c3TkyBFL386fPy/bnFq1atGIESOoSpUqlCtXLtnlkBEEjEAAQpkRvIg+gAAIgAAIgIA+CUAo06ffYDUIgID2CUAo89BHXBSaMGGCrFIXLlygokWLysqrlUx//fWXbHGvUKFC9Oeff1K2bNm0Yj7sAAG/EoBQ5lfcaAwEQAAEQAAEQMCGAIQyGxg4BAEQAAEFCUAo8xAmn8Y4bNgwWaVOnTpFZcqUkZVXK5l4ZFnJkiVlmcM3NeA7aiKBgFkJQCgzq+fRbxAAARAAARAIPAEIZYH3ASwAARAwJgEIZR769ZtvvqFu3brJKrV3716qXbu2rLxaybR27Vpq06aNpDnh4eGWKZuZMmWSzIsMIGBUAhDKjOpZ9AsEQAAEQAAEtE8AQpn2fQQLQQAE9EkAQpmHfjt+/DhVqFBBVqlVq1ZRu3btZOXVSiYeLcej5qQS36SAb1aABAJmJgChzMzeR99BAARAAARAILAEIJQFlj9aBwEQMC4BCGUe+jYpKYmyZMlCcnaLnDFjBvXv39/DFgKbnUfA7d+/X9KITz75hEaNGiWZDxlAwMgEIJQZ2bvoGwiAAAiAAAhomwCEMm37B9aBAAjolwCEMi98V7lyZTp69KhkSS4kcUFJL4mLf3xhfjkiIBfTatasqZeuwU4QUIUAhDJVsKJSEAABEAABEAABGQQglMmAhCwgAAIg4AUBCGVeQOvduzfNmzdPsmSvXr1k5ZOsyE8Z9u3bR3Xq1JFsja9L9vDhQwoNDZXMq8cMx44do8GDBxN/fPLkCaWmpuqxG6rYHBYWRpGRkVSkSBHL7q+vvPKKKu3opVIIZXrxFOwEARAAARAAAeMRgFBmPJ+iRyAAAtogAKHMCz9s2LCBWrVqJVmybNmy9Ndff0nm00oGHv320UcfSZrTuHFj2rp1q2Q+vWZ4+eWX6c8//9Sr+X6zm++OeubMGb+1p8WGIJRp0SuwCQRAAARAAATMQQBCmTn8jF6CAAj4nwCEMi+YJyYmUr58+SxRVVLFL1y4QEWLFpXKpon7devWJb5Tp1SaMGECDR8+XCqbLu8nJycTj5oSBEGX9vvb6AcPHlD27Nn93axm2oNQphlXwBAQAAEQAAEQMB0BCGWmczk6DAIg4CcCEMq8BN2jRw9atGiRZOnp06fTgAEDJPMFOkNsbCzlyJGDuAgolX799Vcy8pQ7Hgl46tQpKQymv1+wYEG6evWqqTlAKDO1+9F5EAABEAABEAgoAQhlAcWPxkEABAxMAEKZl87dvn078SmIUkkv0xT5VMqmTZtKdYcyZ85MPIooJCREMq9eM+zcuZMGDhxI58+fl7WxgV776YvdPEqSi8AtW7b0pRrdl4VQpnsXogMgAAIgAAIgoFsCEMp06zoYDgIgoHECEMq8dBCfmsfXsjpx4oTbGvg0vjt37lDWrFnd5gv0za5du9KSJUskzeBi2ubNmyXzGSUDFwXlRNkZpb9S/QgKCiK+mQMXTJGIIJThWQACIAACIAACIBAoAhDKAkUe7YIACBidAIQyHzy8bt06at26tWQNI0eOpLFjx0rmC1SG+/fvE59GFx8fL2mCkdcnk+w8MoCAAwEIZQ5AcAoCIAACIAACIOA3AhDK/IYaDYEACJiMAIQyHx1euXJlOnr0qNtaIiIiLLsDPvvss27zBerm1KlT6X//+5+s5g8ePEjVqlWTlReZQMDoBCCUKevh69evW6J0+eOjR48sG6bwx+joaAoNDbWso5gzZ07KnTu3JaL3xRdftFxX0orHjx9bftjgUcOjR4+mLFmyKFk96gIBEDAQAYxZBnKmTrsCoUynjoPZIAACmicAocxHF8ld26tLly60ePFiH1tTvjj/Mvj888/T6dOnJSvnU+4ePnyo+BdTyYaRAQQ0SgBCmfeOefLkCe3YsYP27NlDf/75p+X/3r17HlUYHh5OFStWpDZt2tAbb7xBkZGRHpV3zMwFOb6u5KFDhyy3Zs2aRX379nXMhnMQAAETEsCYZUKn66DLEMp04CSYCAIgoEsCEMoUcFvPnj1p4cKFbmviazvxL19VqlRxm8/fNydNmkTvv/++rGYbNmxIfBMDJBAAgTQCEMo8eyZcunSJfvzxR9q0aRPt3btX0c0ygoODqXnz5vTxxx9bxDPPLCOLLXwNRi7cWdP3339PHTt2tJ7iEQRAwGQEMGaZzOE67C6EMh06DSaDAAjoggCEMgXcFBMTQxUqVKCzZ8+6ra1QoUJ04MAB0soUTP5FtUGDBpSSkuLWbutNvs4aX28NCQRAII0AhDLpZ0JycjJt2LCB5syZYxHaeRSr2ql9+/bEo8H4FE05KTU1lTp06ECrV6+2y86nVfkapWZXIU5AAAQ0TwBjluZdBANtCEAos4GBQxAAARBQkACEMoVg/vbbb1SjRg1KSkpyW2PZsmVp//79xNfZCWTiXwD5lKVbt27JNmPfvn1Uq1Yt2fmREQSMTgBCmWsP3717l2bMmEHz588nPt54kngEbpkyZahIkSKUK1cuyz/fOZivH3bjxg3L/4ULF+jq1asuq+U/SKxcuVLWmop8eiUX8mxTqVKlZE1Jty2DYxAAAf0SwJilX9+Z2XIIZWb2PvoOAiCgKgH26z6SQgSWLFkisOk/PFzC7T9bDF9gUWgKtep5NefOnRPKlSvn1kbHPrC1gISEhATPG0MJEDAwATaF0OXriO0Qa+Ceu+4aW3xfYIvgC2wRfJdsHMcXPm7Wq1dPYFGrwrZt2wS2FqLrBmzuHDt2TGBRrgL7AUK0LbYBgPDll1/alHA+/Oijj0TL9urVyzkzroAACBiOAMYsw7nUVB2aPn266HsYf59lPxaZigU6CwIgAAJKEiAlK0NdgjBv3jyBRUO4fNOyfkFk0VzCqVOn/I5s48aNQvbs2Z3sK1++vFu7+ZdYJBAAAXsCEMqe8oiNjRUmTpwosAgwp/HFOu45PrI1GwW2667AIs6eVuTlEfdF8eLFRdvu16+faK1seqZofm7n0qVLRcvgIgiAgDEIYMwyhh/N3gsIZWZ/BqD/IAACahGAUKYC2WnTprn88mX7RTEiIkJgU5NUsMC5SrabmzB8+HBRMYwt0i/ExcUJbFMCl3aPGTPGuVJcAQGTE4BQlvYEYLtXCsWKFXM5ftiOezx6jK0hJvBoMKUTH8f4WJUxY0YnWz799FO75latWuU2AvjKlSt2+XECAiBgHAIYs4zjS7P3BEKZ2Z8B6D8IgIBaBCCUqUR2xYoVQubMmZ2+rNl+YbQeN27cWGA7YqpiCdvOXOBTwFxFeXTq1Cl9SiWfksmnKlntsn1k66qpYh8qBQE9EzC7UHb//n2he/fuomOG7fjBj7lA1rVrV79E0m7dulXgP0Q42rB48WLL023Xrl2iYpo1P49MQwIBEDAeAYxZxvOp2XsEoczszwD0HwRAQC0CEMrUIsvq5VMrX3jhBacva9YvY46PfO0yLrCxDQF8sort4Caw3TWF9957T8iTJ49o+2FhYcK4ceMEntc28agLR7u44JeYmGibDccgAAKMgJmFMj6NO3/+/E7jheP4wc95tBnbDMSvz5mdO3cKmTJlsrOPi2fLli0T2MYAdtcdbebiHxIIgICxCGDMMpY/0Zs0AhDK8EwAARAAAXUIYNdL9g1JzcQW7afBgwfTggULKCUlRVZT7MsnsTXBLDtM1q5dm/hOmXwXOFeJ1/vXX38R33nz8OHDxL6807Vr11xlp8qVK9PXX39NbEF/pzy8Lr57J4twS7/XqlUrWr9+ffo5DkAABNIImHXXy8mTJ9OwYcOICe2ST4V33nmHvvjiC2KCu2RepTOwjQGoadOm/Achj6r+9ttviUW/eVQGmUEABLRLAGOWdn0Dy3wjgF0vfeOH0iAAAiDgigCEMldkFL5+8uRJGjp0KG3ZssXjmln0A0VGRhKbPkm5c+cmtpscsR0o6c6dO8TW0aGrV68Si/iSrPfFF1+0fLmNioqiDBkyuMx//vx5i1h269YtS565c+cS/7KLBAIgYE/AbEIZi3alvn370sKFC+1BiJwVKFDAko8LVYFMbCF/mj17tkcmXLp0iQoXLuxRGWQGARDQHgGMWdrzCSxSlgCEMmV5ojYQAAEQsBKAUGYl4afH7du3E9vljdhCssQ/wKmd2LpAlsi0IUOGUIsWLdxGptnawoW9Pn36EFsc22Jrjhw5bG/jGARAgBEwk1D24MEDatu2Le3evVvS92ztQ5o5cyZpYdx4/PgxsSnwlh8VJA1nGbhAxoUyJBAAAX0TwJilb//BenkEIJTJ44RcIAACIOApAQhlnhJTKP/Dhw8t0xnZzmvExTM5EWFym86WLRs1adKEWrZsSc2aNbNEockti3wgAALyCZhFKDt79qxlPDlz5oxbOGy3Scu0bi6UaSmtW7eOWrduLcukt956i7755htZeZEJBEBAmwQwZmnTL7BKeQIQypRnihpBAARAgBOAUKaB5wEXzX788UfLP5/2eP36dcu0Sjnr6rBF+YmvacanOVWvXt3yZbZmzZrEdq/UQM9gAggYm4AZhLI9e/ZYIsnYbnFuncmnhPNxrG7dum7zBeImX0utRIkSdPHiRcnmFy1aRN26dZPMhwwgAALaJIAxS5t+gVXqEIBQpg5X1AoCIAACEMo0+hzg0zJv3LhhEc24cMb/+bpkbBdLy3plfM0y/p8zZ07Z0yk12lWYBQK6JWB0oYyLRnwKttQ0cT4ubd68mSpVqqRZX44fP54+/PBDSfv4jxVsl07JfMgAAiCgPQIYs7TnE1ikLgEIZeryRe0gAALmJQChzLy+R89BAAR8JGBkoYzvjNu9e3dJQnyq9/79+4lvFqLldPPmTXruuecoOTnZpZn8/uXLl13exw0QAAHtEsCYpV3fwDL1CEAoU48tagYBEDA3AQhl5vY/eg8CIOADAaMKZXv37qVGjRpJRpLxqd98J9969er5QNF/RXmf+EYqrtKbb75JS5cudXUb10EABDRKAGOWRh0Ds1QnAKFMdcRoAARAwKQEIJSZ1PHoNgiAgO8EjCiU8amHr7zyCt27d88toKCgIPruu+9Iawv3uzN69OjRNHbsWJdZ5s+fTz179nR5HzdAAAS0RwBjFsYs7T0r/WcRhDLfWUdHC3TjVio9fChQlixBlCd3EOXKGUzBwb7VffVaKvF/nqpUCqEMGXyrD6XdE+B+vH4zlR49EqhggWDLP/uoigQCXhOAUOY1OhQEARAwOwGjCWV8Y5Fq1arRP//8I+nagQMH0rRp0yTzaSnDpk2bLBueuLKJ7+pZsmRJV7dxHQRAQGMEMGZhzNLYU9Lv5kAok4/81u1U2rojmTbtTqG//hXor4esbAz7TxGpgwssGYnyZCYqkp2oRMEgqlstmOrWDKFSJaQVr/h4gXLWi6O4B2l139wXTvny+qi8iZhp1ktsjyY6+GsS/bAxhZbtS6U7nLOjH5mbInMSlX8uiMoWD6KuHUPopXIhZkWGfntBAEKZF9BQBARAAAQ4ASMJZXztrubNm9P27dslnVu2bFk6evQohYeHS+bVUoY7d+5Q3rx5RU3iOwdfu3ZN9B4uggAIaI8AxiyMWdp7VvrfIghl7pkf+CWZftyaTKv2ptLFG+7zyr4bQdS6fBDVZ8JZvVoh9MLzzsLZ0NHxNHllWjQZr/f8loxUrIhzvtw1YukeF+xUSgP+E0TTP2MG+5iOHkumyp0TfawlrXibqkG0+mvvbOIC5OQZCTRqKWMb77k5dV8gGj0wlOrVCfW8MEqYjgCEMtO5HB0GARBQioCRhLJ+/frR7NmzJdHwdckOHTpE5cuXl8yrxQx8wf6rV686mRYVFUXLly93uo4LIAAC2iSAMQtjljafmf61CkKZOO8ff0qkD75Mpr/9sT8PizrrXTuIqpQPJj7Vb/WWVPrpqGBn2Ik1GalcWWehLKhqLFG0XVZFT3o3D6I5U7wTpWwN+e1oMlV9Uxmh7KUiRH9syWRbveSxwHB+tyKRukxmGzI9kcwumaFUQaIFn4RSrRoQzCRhmTgDhDITOx9dBwEQ8I2AUYQydx+0HQl99tlnNGzYMMfLujmvUKECHT9+3MleLhL26dPH6TougAAIaI8Axiyy/LCBMUt7z01/W+TutbBy5Upq3769v00KWHt8Ot6qtYn0/lfJdPlmwMwQbfjQ0jCqWtl52t/S7xPptz9S6a/zqbTzJCv6NAhNtB63F9nMzldLEL1QNIgK5g+inNmDqFqVDKLtuq1H5Ob9+6m0dGUSXb4m0CX2v/oPplzxaatyE7OtIdscvVyJYKr9ajC1/k+Y3JKUkCBQpz7xtPYXe/FRdgWuMjJRc8EHGahHVzbHFgkERAhAKBOBgksgAAIgIIeAEYQyPtWyWbNmlJLiuLiDM4HatWvT7t272QK3+l1no379+pY+OPbu1KlTVKZMGcfLOAcBENAYAYxZaQ7BmKWxJ2aAzIFQlgb+jxPJ1KxfIt2444Ej2OoRveoFUeWXg+lZvvh7JBOY2GOWzEF0645AN9nC8FevMwFrfwrN2sdEmkce1O2Qdfe8UKpb23300j9nUqjzoAQ6et6hsIzTRmwq6NKvMlLePP75fMYjvPb9nER1hyRJcnm1FNE3X2Sk0iWdI+qkunb/gUCN3oyj3y9I5fT+/nttg+mLseE+b97gvQUoqVUCEMq06hnYBQIgoHkCehfKHj9+THy9MbGpiI7wIyIiiH8xK1y4sOMtXZ23bduW1qxZY2dzvnz52Adijf38bGchTkAABDgBjFlpzwOMWXg9WAlAKCP6am4CDZzGfuyTE43FtKq36wZR57YhVLdWqEc7UZ46nUKLvkuiyatZQ9K/LVpdZHncOC2UWjRxL5TxjFyA6j4wjr7Zzg48SI+PRFBmJvD5O/17OZWKNGWLhblg3/Blom3fZ7JMSfXUNr4LaamWcXTnrqclPc8/4s1g+nQUU06RQMCGAIQyGxg4BAEQAAFPCOhdKOvfvz/NmjVLVpf5dEs+7VLvqUePHrRo0SK7bvCpKXyKChIIgIC2CWDMSvMPxixtP0/9aZ2ZhTIebdS5fzxt/l2eqNSpThBN+zQj5cntW9TV9Rup9NHnCbRgK2tXXtO0YnwIdWgjb7oh21uJXmoWS6euyHwmse4If3u25pfMmmVlK1o3li6J/dbIuntndzjlzuU5bz7JoV77WNr/tywTFMm0Z34o1WHiKRIIWAlAKLOSwCMIgAAIeEhAz0LZgQMHqFatWuzXS+lPedmzZ6cLFy5Qjhw5PCSkvexiX7RnzJhB/DoSCICAdglgzHrqG4xZT1mY/cisQtmfJ5Pp5bfY4vIy1smKYB9dNk4Ipfp1lRVBjv2RTBV7MBtkLC7/9cgM9HZn+Wthnfw7hV5sl+AyUsvxeX9keRhVquC8BppjPjXOC9SKFZ3yOmVABhrUX36fbW0bPDKepvwgHqZWqThR41eCqXTxICpZPAMVKRxMSWwG6O07qXTyVApt2JlKa/lmCp7uPfAM0b0dEZQzh/8j82z7jmPtEIBQph1fwBIQAAGdEdCrUJaYmGjZtZJPpZSTxo8fTx988IGcrJrP06tXL1qwYIGdnSdOnKBy5crZXcMJCICAdghgzMKYpZ1no7YsMaNQduFSChVvzUSkOGlfNHiJaNOSCMqYUR3x49btVGryVjz9cdG9LTMGZ6D+vTwTjQYMi6MZ66V/zOQtD+8UTBM+8v/UQb6BQoZybOdOR02LbbQZ90sEhYd7zt0igrZnKpdD10s/y3a+nCRPEIyPF2jm/EQaspCFprGZoXLTu22C6cvx/uco1z7k8y8BCGX+5Y3WQAAEDERAr0IZj0YYMGCALE/kz5+fzp8/T5kyBS6sX5ahMjNFRUXRihUr0nPnzp2bbt++zdbP8PzDXHolOAABEFCVAMYsjFmqPsF0XLnZhDIuTBX8TzylyFhYvz7bZXHzdxEUFqbu+3tiokAdesXT+kMOyo7N82pi32Aa+q5nAgzfaTJXHabysGgpyZSZLZt2OJPfF6TnmyiU56KWQ/q0VzCNGOxZf61V1Hgtlg6etp6lPfZuHkRffR5BoR4GBUZHC9SmZxzt/NO+PpdnTMtM+E3954zL9nFDUwQglGnKHTAGBEBATwT0KJTFxsZS8eLFZS9eP3PmTOrXr5+e3OLW1gYNGtCuXbvS87Rp04ZWr16dfo4DEAABbRHAmIUxS1vPSG1ZYyahjIseL7eKE18Py8Et9ViQOBfJ1Iokc2iOuFhWq20cHT7reCft/KO3g2nMcM+FozGfxdPH3ziGa4m3seHLUGrZ1EMlSbwq2VdHfBJPE5Y52Mc2t7y3P5xy5vR8bbLde5Oofm97ZfDtRkH09VcsRM3LxKPeeHTerA2uhUzbqpePC6GodvLWk7Mth2PjEYBQZjyfokcgAAJ+IqBHoWzixInEF+aXkyIjI+nff/9lv+D594OXHNu8zVOqVCk6e/bpJ9lp06bRwIEDva0O5UAABFQmgDELY5bKTzFdV28WoYyLHdVfj6VDZ6Td9Ww+onPb/CeSWS169EigZ5vG0ZMH1itPH4d0CKZJn3gulHFxMFsNNsfUXjt6WrHNEd9hcvsK/0b/Z6seS9H3bYxgh281DKJvZngnbHVg0V+rfn4qaPH1yA5v8D1Sjj9/GkfFyoosq1WWaN8a/3K0J4gzrRCAUKYVT8AOEAAB3RHQm1AWHR1NxYoVo3v37sli/eGHH9K4ceNk5dVLpoIFC9L169fTzT19+jRx8QwJBEBAewQwZhFhzNLe81JLFplFKJuzKIH6TmTrTUklNsvy8NIwqlIpMAvb7z+QRLV7OKtafVsG0azJ3olHo8fH09jFDlFbLjic/SkjlSjGQrr8kCxribVznnbp7cYCXGjMXp2JglY3Mxde2JSRihZWpj8W0bEOq5/9uU3MTcIxCGVuGZnkJoQykzga3QQBEFCegN6Ess8//5yGDx8uC0RwcLBlbbIiRYrIyq+XTPxLxZQpUyguLo5atGhBCxcu1IvpsBMETEcAYxYRxizTPe096rAZhLL7DwTKVZepG2z9fqn039fS1rKSyqfm/Z7vxtHCrU+jonhbXRsE0bczvRPKLGuV1WJrlVkFJDfG9/tPEM2c5F07bqoVvfXeiHiatsZewCsaycSt3d6JTPO+SaDenz3tpDfruokaanNxyowEGjzjaRs2t+wO4457txGBXSU40T0BCGW6dyE6AAIgECgCehLKBEGgEiVK0IULF2Thatq0KW3evFlWXmQCARAAAaUJYMxSmijqMyIBMwhlvf8XR/N+sheeRH2ZhejJngh65hl1F+8Xbdvm4s1bqRRZjwlbNhpS2+pB9MMi7wWsdwbF0fzNMhiwKKyHByIoWzZ1GSQkCBRe1Vm89GZ3Tyu6/3SJo42//X8f2eYE8awfSq8xx3fDjOB2OwfCWc2wPP69PiM9X1qZSDa7inGiKwIQynTlLhgLAiCgJQJ6Esq2bdtGTZo0kY1vzZo11Lp1a9n5kREEQAAElCSAMUtJmqjLqASMLpQd/zOZKnRkqoYMjWjRyAzUrTPbtlAD6c0+cbRsz1Ojm1YIos3LvRfKzpxLodItZYTUsb5P6B1Mwwd5vh6aJ9iWrUykN0cn2xdhutKjXyIoa1bPRboUFuQVUj42fS22BSMyUI+u6viyXfc4Wn3wqW/sO5F2tnlGKDVtaJz1ecX6iGvSBCCUSTNCDhAAARAQJaAnoaxt27bExS85iS/if/nyZQoJCcwaH3JsRB4QAAFjE8CYZWz/onfKEDC6UFardSz9fEoGK6ZpxB6OoIgIz0UaGbV7nMUi8HV4GrZU83mi/Wu9m5JobbxZpzjacsy9wGPJ+wzTmw5lYp/hrCWVf6zaMpZ+O2dfb+d6QbRktndi4K+Hk6la1//nxTAlMvvV2kfqm6UJ1G2c++mXW2eFUuP6EMrsPWy+Mwhl5vM5egwCIKAQAb0IZTdu3KBChQpRcrLDr38uOAwYMICmT5/u4i4ugwAIgIC6BDBmqcsXtRuHgJGFsitXU6lQQzaFUUbq3TyI5kzxTqSRUb1XWZ6tHUvXbqcVfakI0R9bfBPK9v6cRHV7Om8UIGbckjEh1DkqTOyWz9f+OZNCz7dyjm47+G0YVXvFO3Vu/Bfx9OH8tLmqH3cPptHvqxcR9/vxZKoU9VTEFANyhm0iULI4pl6KsTHTNQhlZvI2+goCIKAoAb0IZXPmzKG+ffvK7vvWrVupcePGsvMjIwiAAAgoSQBjlpI0UZeRCRhZKPt8ajwNn2uz0JcbR/7xQxi9VM47kcZNtT7dst2tskh+oot7fBPKuDGlG8TSmWvSZpUoQHR2l+/tibX07gfxNH2tvV/y5SG6ud/79k6fTaFV65Po6EmB5n+RkXLnChZrWpFrV6+l0nMN3AuwiSfUi2hTpBOoxC8EIJT5BTMaAQEQMCIBvQhlrVq1og0bNshyQebMmenevXsUFqbOL5GyjEAmEAABUxPAmGVq96PzHhAwslBWuG4sXb4pDaNUQaLTO70XaaRb8C7HX6dSqFzrtMirzDmIHv/iu42ia4O5MG/3vFCqW1vZ6YOuFvGfMiADDeqvzppiLrrn9eVbt1Mpf203Qhmbuioc9d1XXhuIgpohAKFMM66AISAAAnojoAehLCEhgXLlykUxMTGy8LZp04ZWr14tKy8ygQAIgIDSBDBmKU0U9RmZgFGFMluRScp/g9sH0+Sx6k3Vk2rf3X2+W+X9R0QvPx9Eo4b6biNfQSP0Vbbo/RN3rabda1Q+iLZ9r+x0VFGhjgV/3WM7VObMoY314aTIHD6STK90dj31skJRot83QyiT4miG+xDKzOBl9BEEQEAVAnoQyjzdOW7RokXUrVs3VXihUhAAARCQIoAxS4oQ7oPAUwJGFcqGjYmnid/bT+972mv7ox8+D6G2r5knCn7StHh6f7Y8Nqc3ZqRSJZRba0tsEf/2NYNo5QJlBTl7Dyt7Jir22TTRrXEQLZqun/7YmI5DhQlAKFMYKKoDARAwDwE9CGWDBw+mKVOmyHbK9evXie96iQQCIAACgSCAMSsQ1NGmXgkYVSgr3yyW/rgozytXd4VTwQLqrWklzwr/5YqOFihb9TgiGfsz9WoWRPOmKiP6uFrEf8/8UKpTS9kpnmrSHPNZPH38jWuhcc8C1p+a+umPmqzMXjeEMrM/A9B/EAABrwnoQShr0KAB7dq1S1YfCxcuTJcuXZKVF5lAAARAQA0CGLPUoIo6jUrAqEJZ0EtseqHr2XFP3WnS9aQGDo+jr9YJTzm4OmLBZPf2h1POnL4LiWKL+Cu19por89W4/mafOFq2R5ydHvujBiPUmUYAQhmeCSAAAiDgJQE9CGV58+alO3fuyOphx44d6fvvv5eVF5k8I8AFyL1799LVq1fp1q1bdPv2bYqLi6OsWbNS9uzZKUeOHPT8889TxYoVqWTJkhQc7PuHWs8slJ87KSmJPvzwQzpz5oxl/buRI0dS0aJsUQ8kEFCAAMYsBSAqUIUexyw92uyrq4wolN28lUqRddwstm4DrUXlINq4VJmIKZtqNX946XIKFW3CNgoQ13vs7P+4ezCNft+39dFcLeL/eZ9gev893+q2M9YPJ+52Dp3QO5iGD9JXf/yAzLRNQCgzrevRcRAAAV8JaF0o44JM/vxsT3KZiU/RHDRokMzcyCZF4Ny5czRjxgzavHmzRVSSym+9ny1bNurcuTP997//pTJlylgva+KR74jarl072rNnT7o9o0aNok8++ST9HAcg4C0BjFneklOmnB7HLD3arIy30moxolC2a08SNeiTJAvTwNbBNG2COYWN1m/H0bpfZShlbF36hF8j2G7m3i+2/92KROr8kcNcT/Z73q294ZQ3j3Z/2HN8Ep05l0KlW6btROp4jxieG3vCKX8+/fTHqQ+4oCwBAQkEQAAEQMArAps2beKfUET/J0yY4FWdShbavn27qG2ubD548KCSzZu2rtOnTwtM6BIyZMjgEX9HvwQFBQmNGzcWNmzYIKSkpASUZ2pqqrBgwQKBRfs49endd98NqG1o3DgEMGYFxpd6HLP0aLMa3p0+fbrTmGx9L1m5cqUaTape57TZ8QKVjpH1P+azONXt0WoDvx5OksWIs/x6SbxP3ajSwtkfr3WN9anOQBQeOjrOJbPmb+ivP4FgaKY2EVHG3k2QQAAEQMAbAlqPKJs1axb1799fVteYKEMxMTEUEWG+KQyyAMnMNH/+fGLCkWVapViRzJkzU6FChShjxoyWKZg8goaJYGJZ7a4VK1aMZs+eTUw4s7uu9snly5ct03GXLFlCJ0+eFG1uwIABxL6sid7DRRDwhADGLE9oKZNXj2OWHm1WxlvOtRgxouydQXE0fzPX/6TT9EEZaEDvjNIZDZpD7qYHhdjkgn/3sNAyL5KrRfy3zgqlxvX1s+g9/6gVUoWtfcf+nBILIju1LiOVKaXcDqFObeCC/giYSRVEX0EABEBASQJajygbN26cy1+a2buV3T22kL+SaExXV2JiovDmm2/aMbUyrlq1qsAEAIGJYk5cYmNjhbVr11oi0Jh4JlreWg9bt0xQO1KRrZ0m8MhC9uVLqFGjhsCj2qztu3pkQplTv3ABBLwhgDHLG2reldHjmKVHm73zjvxSRowoa9Yp1mXUj2Ok2dLvE+TDMmDO1esTZLPavivRKwJs4wDnNqrECCzQXFdp/SbXrD6aYN7IRF050c/GIqKMffJHAgEQAAFvCGg9ouz999+nSZMmyepao0aNaNu2bbLyqpnp7t27dOXKFWLvhcSmLlJISAiFhoZaIrDCwsIsjzwaix/z63ITmzpoifLiC+jz//j4eGJfuogvTM9TqVJvEgf6AAAq3klEQVSlKFMm735t5bZ26dKFvvvuOztzXn31VeKRD+XKlbO77urk/Pnzlmg0JsC6ymK53r59e1q2bJmFjduMEjf5Ivy///67hUV0dDTxtX4ePXokUcr5NiLKnJngincEMGY95YYx6ykLfqTXcda+F8qfGTGirHabWNr/tzxWm6aHUvPG8j8LyKtVP7nYRxsKqxZLKTLeuuuxjyK7fvDsc46rRfw/6RFMo4bqa2245m/E0ebf+W9/9ilfHqIruzKxz5T213EGAvyNBwkEQAAEQMALAlqPKHvnnXcko4HY26AlD5ui6QUBZYv88MMPAo+astok9cijncLDwwW2+L2QJ08eoUCBAgKPjCtSpIjw7LPPCvny5bPcY6KaZJ0lSpRgv4569/PokCFDnOrn63bx6Adv0vjx453qc2QxdOhQb6q2K1OnTh3JdqztuossQ0SZHVac+EAAYxbGLFdPH72Os676o9R1I0aUVWjqvB6WYySZ9XzXHu/eZ5Xir4V6PFnT7a9TyR6ZzCP2rKzTH8vECNdvBHbdVI86wTKfu5AsELM7vQ82a+D9cijJ0+qQ3yQEIJSZxNHoJgiAgPIEtC6UdezYUbYQwna8VB6QhzW+9tprsu21CjhKPj548MBDiwXLNEVbEYkff//99x7X41hg8uTJkiz4lE1fkpRQxoXHvn37Cr/99puwaNEil/ZAKPPFCyhrSwBjlv2UeKnxzSxjFp8Ortdx1vb5rcaxEYWy0g3EBQ0xkWPtBnNPveTPqSdPUgV6UR6zt/p7tmC92CL+TaM8q0ON572nddZrK86n72D99cXTviO/9wQw9ZJ9CkECARAAAW8IaH3qZatWrYjtmCira3wqX6dOnWTlVSsTE25o3759alXvtt4cOXLQ/fv33eZxvMmnRlWuXJmOHTuWfosJS8QXJPc1sbd1Ymub0ZEjR1xWxSLp6OzZs8Si6VzmcXdj0KBBdOLECcsU1uzZsxPb0ZLy589PL774IlWqVIkiIyPTi3/99dfUvXv39HPbA0y9tKWBY18IYMyST88sY5bex1n5HvUupxGnXr7aKpYOnZHHY8GIDNSjq3kX87dSGjYmniZ+z+ZhSiW2Vv2tPeGUNw9bvV4iuVrEf+O0UGrRRD/zFFevT6R2w5KdehuRg+j29gjKnDnI6R4ugICFgPcaG0qCAAiAgLkJaD2ijAlfLqOA2BuA3b1du3YF3Jnr1q0TqlevLpQsWdIyffK5554TmFgjMBFHkFro3rE/UudsfTMhV65clumaFStWFJYvX+5x/7m9tu3waZ+PHz/2uB5XBZhoaFe/bVvWY7UX97fahogyKwk8qkkAY5b9uGx9nfNHs45ZZhpnvXltGTGi7PW35C/m//lULMLOnzdXr6UI9Lx41JRjJN6oT+UxE13Ev1KMkKKjWZeWaLuKIlxeiBH+/sezaajevD5RRt8EMPVS3/6D9SAAAgEkoHWhjE+bs/2i5e6YRRYFkKS8ptki/MLVq1eFPXv2CHxqYvHixWX3j/e9Vq1awubNmwW2YYC8BiVyOU4Tq1evnsCiDAUWmSZRUv7tdu3aue0jF+dS/PCpFUKZfJ8hp/cEMGbZC2UYswTBTOOsN68cIwplvQfJF8qGjpYn+njDVm9lOr0jk9vLMUJ8vPs1Wfl9eslZYBo5Tl+8//ehyI6dbH2ydRsxZVdvz+9A2AuhLBDU0SYIgIAhCGhdKPvggw/ciiy2wtmFCxd055OTJ0/K7l+hQoWEhw8fKtZHtmum8Mwzz4i2z6YwCrt371akLTa1UrQNW99t3LhRkbbcVQKhzB0d3FOKAMasp0IZxiyBfZk31zjrzevIiEIZF2Mco6BcnXcbgDWmrM+bY38kyeY29+t4azHRR9FF/JnA9O8V/YSTnfxbfAH/8V/oS+wTdRAu+oUAhDK/YEYjIAACRiSgdaHs888/lxRZrILLqVOndOciOSKStX8rVqxQtH+nT592yzZLliyCNwttixnJhTdrP8QeR48eLVZM0WsQyhTFicpcEMCY9VQow5glCGYbZ128LNxeNqJQNn1OvGzBp2VnCGW2T5DqrZyjwMRExnw1Y2yLOR2LLeJfv537Mk6VBPBCQkKq8EIjZxb9h+L5EkC36K5pCGW6cxkMBgEQ0AoBrQtlbAF2twKLreiyd+9erWCVbceaNWtk9S88PFyIiVH2A96OHTsk254zZ47svrjLWLduXbdtdenSxV1xRe5BKFMEIyqRIIAxK00ow5iV9kQx2zgr8fIQvW1EoWz5qgTZQlnxesq+t4tC1tHFTVsTZbP7aVuiaM9OnWaRWCx6zPF/zY/6ma7YuY/zNNS23WKFVPczTkV54KJ5CUAoM6/v0XMQAAEfCWhdKGO7MboVWGyFslWrVvlIw//F33//fVn942uHKZ2+++47ybaHDRumSLP9+vVz21bNmjUVacddJRDK3NHBPaUIYMxKE8owZqU9o8w2znrzOjKiUGaZMici1DgKN5bzMjFCdDTUD+tzhwtBuao7i1xi7Gq+Li4yii7iXz5GSEqytqLtx0nTnCMS67aJERLFdUFtdwbWBZQAhLKA4kfjIAACeiagdaEskX0qkLtb5IwZM3TnCr7Qta3Y5+r4k08+UbxvfNF+V+1Zr48bN06Rdt19EeJt8bWM1E4QytQmjPo5AYxZaUIZxqy014PZxllvRgF37w8rV670pkpNlCG24LyYuCN2bfN2KCC2TuPrj4lxErv250l79cuyiP+Lzuzf/0gf63pt3MIi6ph4atvXqi1jhLg4iKm2zxEcyyMAoUweJ+QCARAAAScCWhfKuMGVKlWSFHS42DJq1Cin/mn5Al/kWa4IuH//fsW78u+//wpBQUFu2W7ZskWRdhcuXOi2nRIlSijSjrtKIJS5o4N7ShLAmEUCxqy0Z5TZxllvXkdGFcpe6+o8dc5W/LA9HjZGHyKON/71pgwXheQKjW+yKYq2ydUi/ucvJttm0+TxX6fYlNFy9iJZWbZO2ePHEMk06TAdGAWhTAdOgokgAALaJKAHoeydd95xK7JYo5969eqlTcgurOJrqlltd/eYKVMmS5SKi2p8uty5c2eXNrzyyiuKtfvZZ5+5bIf3vVmzZj71Q05hCGVyKCGPEgQwZmHMsn0emWmcte233GOjCmVfzZUfFVWxmfgUQrkMjZhv1Kcydw59Pka4cfPpTpaVmtsLTVyQrPGa9vneu5ciZKhqb3uROjFsUyWIZEZ8fvurTxDK/EUa7YAACBiOgB6Esh9//NGtyGIVmcqWLasr/3z88cey+tW4cWPV+pWcnCwMGjTIElnGF9+2smzevLlw9epVxdrt1q1bet3WNmwf3333XcXaclURhDJXZHBdaQIYszBm2T6nzDTO2vZb7rFRhTJLdJDcdcqY2MOnDCI9JXDrdopAZe2FI9soPNvjDz5Oi8j7+x/xRfy//0Hbi/jzCDpHgS8HW6eNM0ACAV8IBPHC7AM3EgiAAAiAgIcE2Pop1KJFC9FSEyZMoOHDh4ve8+dFtuYP5cuXjx4+fCjZ7IULF6ho0aKS+bSQge0ESSyqTNIUf/jh+vXrdP78eQoODrawZlMhJe3yJEOxYsXo4sWLLovMnDmT2IL/Lu8rcYPtRkjdu3cXrWrAgAHEvqyJ3sNFEPCUAMYs9d879Dhm6dFmT5/73uT/6quvaODAgaJF2Rpl1L59e9F7ergYVDGWiP3JSZumh1LzxqFyspomT7cBcfTNdhlf8zMyzL9G0LCP4+mrdQ752b3EI5koVKNok5OJmr0ZSzv+sHFrNqKra8OpYIFgm4s4BAHPCUAo85wZSoAACICAhYAehDJuaI8ePYhFBEl6jYsdXPTQeoqNjaUcOXIQ/0ItlX799Vdi0yClsmn2/unTp6lMmTJu7fvll1/o1VdfdZvH15sQynwliPKeEMCYhTHLk+eLr3m1Ms560w8jC2Wjx8fT2MWpsrBUL010YH0mWXnNkunU6RQq+1qCrO5+8d8MNHheClPF7LO/1zaYpn4abn9RI2ep7KnRvkccrfnFRtx7hujCmoxUtHAGjVgJM/RMAEKZnr0H20EABAJKQC9C2fbt24lNQZRkxfNs3bpVMl+gM3AbmzZtKmlG5syZ6cGDBxQSEiKZV6sZ+vTpQ3PnznVpXpUqVejw4cMu7yt1A0KZUiRRjxwCGLMwZsl5niiVRyvjrDf9MbJQduduKuWtE0/E9Bs56bdlYVS5on5fO3L66Gmehh1iaeefnpZ6mv+fDRmpdEltik7dB8bR19tsRDKm5/2zSrv2PqWKI70QgFCmF0/BThAAAc0R0ItQxmfYv/zyy3TixAm3DMPCwujOnTuUNWtWt/kCfbNr1660ZMkSSTO4mLZ582bJfFrNcOvWLSpSpAixHT5dmuivqTUQyly6ADdUIIAxSwWofqhSS2OW3O7q0WbbvhlZKOP9dBJDbDvvcNyyShBtWBLhcNXcp7v2JFGDPkleQahUnOjIJm1G6Q36MJ6+XG0TbRhG9AcTSl8qB6HUK2ejkDgBvkYZEgiAAAiAgOcE9LCYv7VXa9eu5T+7Sf6PHDnSWkSTj/fu3RNsF8531ye2Ppkm+yDXqBEjRrj1F1u7TOALXfsjYTF/f1BGG7YEMGbZ0tDHsZbGLLnE9Gizbd+Mupi/tY+nTosvMG+7GL3t8emz/nlPtNqnh0e++6MtI7nHi5dpcxH/MZ857Oj5Qoxw+EiSKq5IYfsBrN2QgN0zVaGr/UoRUSauH+IqCIAACEgS0EtEmbUjlStXpqNHj1pPRR8jIiLozJkz9Oyzz4reD/TFqVOn0v/+9z9ZZhw8eJCqVasmK6/WMvEpo3wRf1ebMAQFsV/ON2xwuZmE0v1BRJnSRFGfHAIYs+RQ0kYerY1Zcqjo0WbHfhk9ooz3t1mnONpyzGaKnSMEm/N2NYJo1UJtRZWdv5hC9+8LVKVSYKKdlixPpK4fs1XvPUksQivucASFhwd5Ukr1vF/OSqBB023m4rJZofvmhVKtGursNjDxy3gaNieVcuQiun9Am9F1qkM3cwPa1/JgIQiAAAhok4CeIso4wS1btriNUGLvhZb7Xbp00STw1NRUoXTp0rL6kClTJoEt9q/JfsgxKioqym0/hw0bJqcaxfIgokwxlKjIAwIYszyAFeCsWhuz5ODQo82O/TJ6RBnv7+UrKQKVkx8VNX1OvCOmgJ3v3J2YZvvzMcKPPwUmQot/FKIK8vnxiLN+Q2IDxsxVw4uWxNtHxjGmW3ao9znvxF8smpG1wXn8933t8XDFCdeVI0DKVYWaQAAEQMBcBPQmlHHvsN3k3AowXCxj0UoCWyBec86cOHGipO1Wsa9hw4aas1+uQUuXLnXbz9q1a/ttyqXVZghlVhJ49DcBjFn+Ju55e1ocs6R6oUebxfpkBqGM93vF6gR7kYSJFy6nEJaJEbbvUk9AEfOD2LXZC5mww2xJt5OJLoH6/W78Fw7TFd3xY/e4SKSltGot878tS3bMp0SqlR4/ThXy1Hjqu7//0RYPtfqNeu0JYOqlmcMJ0XcQAAGfCOht6iXvbExMDFWoUIHOnj3rtu+FChWiAwcOaGYK5t69e6lBgwaUkmITcu+mB2PHjiW23pqbHNq8dfnyZXrppZfo0aNHogaWKFGC9u3bR5GRkaL31bqIqZdqkUW9UgQwZkkRCux9rY5Z7qjo0WZX/THD1Etr37v0jaOlu/nvSDISm4l3em1GKlXC/zs2prI15t8bEUdfrbO3ddp7GWhgn4wyjFc+y8OHAuWoGUckYwbmC4WITm7TzjTDLTuSqNlAtiGBzdr9330SQm90YPNDVUiPHwtUs30c/XkprfKKxYiO/qQdHip0GVW6ImCvm+EMBEAABEBALgE9RpTxvvFosdDQULdRS+w9QyhbtqzAF88PdLp27ZqQL18+SXu5zdZ/JiYF2myP24+NjRWqVq2a3gdrX6yPfNopZxGIhIiyQFBHm1YCGLOsJLT1qOUxyxUpPdrsqi/8ulkiynhfeZRP2CtPo3zSI7VcREdlfjVGuHDJv5FA9+6nCg3aO9s48cs43oWApr6DY59Gt7lgxpkuXKydqav7D7Cpq2Xtec7/Rj37jh5LEp5vaN/ekuXqRa4F9AmBxiUJIKKMfQNBAgEQAAFvCOgxoszaTzbthN566y1i635ZL4k+8sXwd+zYQWzNL9H7al88f/48vf7663Ty5EnZTbFdMS0RWWFh6vzaKNsQDzKyd2tq3749rV69WrQUEy1p165dxARD0ftqX0REmdqEUb8UAYxZUoT8e1/rY5YYDT3aLNYP22tmiijj/T53IYVKdkogEg+6tkWTdswiy5aNDqFO7dX9PMA/Ss39OoH68YXmmXm26ZMewTRqaLjtpYAc800FSjRnxvGf41wltt9ADFvEP1OmwC/i//vxZKrUJZGIBZNZU9acROPeUi5KkPvtwSOBbt8V6OcTAp3419rS/z+yp038bxGUMWPgeThYhlN/EJCU0pABBEAABEBAlIBeI8qsnZk3b55lPTL2XuMyionfq1ixonDq1ClrMb89bty4UciePbuTbeXLl3drd7169fxmo1INDR061KmfVr+8+uqrwu3bt5Vqyqt6EFHmFTYUUpgAxiyFgfpQndbHLLGu6dFmsX7YXjNTRJm13zdvpQjF69lH/UhFl3XuEys8eJBqrULRx30/JwpF6ojbM+KTwEeS2Xa2xZvuo8p6vquNRev5mmD0kjhTKV8reX/AMG3wsPUhjv1HAIv5+481WgIBEDAYAb0LZdwd06ZNcynQWIUa/hgRESHMmDHDLx6Mjo4Whg8fLiqG8UX64+LihJ49e7q0e8yYMX6xU6lGRowY4bIvnTt3FuLj1ZtmILcPEMrkkkI+tQlgzFKbsHT9ehizHHuhR5sd+yB2bkahjHOIjk4VqrfyUEh5IUYYODxOuHotRQylR9fYJtzCgV+ShDZvuxaePvhYWyIZ7+DPB5PcTr/kUw8DnS7+y0Syih761s1UUl+Es1On/Tt1N9Ds0b49AQhl9jxwBgIgAAKyCRhBKOOdXbFihZA5c2aXYo2tYNa4cWPh0KFDshl5kvHJkyfChAkThFy5cona0qlTJyEhIW2tiHPnzrlcZ23//v2eNBvQvMOGDRPtK995dPz48QG1zbZxCGW2NHAcaAIYswLnAb2MWbaE9Gizrf3ujs0qlHEmCQmpwvAxbDdHtpukR2IIy9++R6zA1566dl2+aJbCsu7emyj0+R8Txyq4bjOyZoxFkHLnt0Dee6GRuO0lWJReoNP1GylChqri9nnkYwWEsyotAs8j0P4we/tYo4x9A0QCARAAAW8I6HmNMsf+/vPPP9SuXTv666+/HG+JnvO1y9577z1q06YNhYSwRS28TOxNmH755RdatWoVfffdd3Tnzh2nmvhaY6NHjyYWEUBMQEq/z4Qk+vDDD9PP+QET/Oj+/fvENiuwu67FEzYNiCZPnuxkWrZs2WjJkiX0n//8x+leoC5gjbJAkUe7rghgzHJFRr3rehqzrBT0aLPVdjmPZlujTIzJ3/+kUOdBCXTsothd6Ws5chE1fymIIvMEUV52nCdXEPusQXTrjkDXbhH7F+gKOz58hdUV76Y+tnTWhJ7BNGRAOPtc5CZfgG+tXJNIHUc4b385a2gG6tsjMLtyciT376dSudfi6Ybzx8CAEFNzZ82AdAiNekwAQpnHyFAABEAABNIIGEko4z2KiYmhwYMH04IFCyglhS1IKyPlz5+f2JpgVKtWLapduzbxRedtxSzHKni9XIz77bffiO1kRywqj9hOjo7Z0s8rV65MXKQpV65c+jXrAa+rRo0axCLcrJeoVatWtH79+vRzrR5wzlOmTHEy74UXXqC1a9dSyZIlne4F8gKEskDSR9uuCGDMckVG+et6G7M4AT3a7KnnIJSlEWO/udH8bxOo9zT22SXOU4q+52/4MtGCSeFUuFCw75WpXAP/eBdSLZYo2qYhJvJF/xpBWbI8/THS5q7qh2xCAVVtHUenuBiphcT0wgS2iH9YWGB4aAEBbCCCUIZnAQiAAAh4ScBoQpkVA99hkv8Kv2XLFusl2Y9Zs2alyMhIYtMnKXfu3OxDVxZi0yUtkWJXrlyhq1evUmIi28VIIr344ovEpstQVFQUZcjAPsG5SHxXTC6W3brFfvZlae7cufTOO++4yK2NyzwSj62z5GQM3/WSTXG0RMU53QzwBQhlAXYAmndLAGOWWzw+39TjmKVHm71xFIQye2pJbIfE1esTadKiZPr9gv09xc9Y4HqHV4Ooe1QINWmg/Sh22/5PmZFAg2c8/UG0S/0gWjwrwjaLX4//+34czfyRr0ShjfRum2D6cnzgdyrVBg3zWgGhzLy+R89BAAR8JGBUocyKZfv27TR16lTasWMHJfFPnyqn4OBgS2TakCFDqEWLFm4j02xN4V+S+/TpQ2yRf4utOXLksL2tmWM+zXTgwIHENkWws4lH4I0dO9ZpGqldpgCfQCgLsAPQvCwCGLNkYZKdSY9jlh5tlu0QkYwQykSg/P+lE38l08yFSbT4F4HiHrjO59EdpodFVQuiTq9noKYNQ3UbccQjuLJUY6F3///R7tclYfRKlcDNFx02Jp4mfp/qkSvUzHx6Y0YqVcL1j7Rqto26tUMAQpl2fAFLQAAEdEbA6EKZ1R0PHz60TGfk64jxL6JyIsKsZaUe+XpcTZo0oZYtW1KzZs0sUWhSZfR4n39569+/P82ePdvOfB5xt3TpUsuUUbsbGjuBUKYxh8ActwQwZrnFI+umHscsPdosyxluMkEocwPH5tbNW6m09+dk2n0ghXYeF+jcQ3aTrzfmSpthmlFEFqIXcxMVyhtEhSKDqG6NYEvkmFGm4w0ZFU9f/JBKhfIR/bsnkw0tHIIACHACEMrwPAABEAABLwmYRSizxcO/gP7444+Wfz7t8fr165ZplfwLilTii/LzNc0KFChA1atXt4hjNWvW1MXC+1J9c3efs+ERb/PmzbPLVrRoUdqwYQPxdck8TVy4GjNmjKUY24mU5s+f72kVHuWHUOYRLmTWEAGMWZ47Q49jlh5t9twzziUglDkz8eRKdLRAt++m0t27AlsmQqC8eYOpYGQwZc2Ktak84Yi8IGBEAhDKjOhV9AkEQMAvBMwolImB5dMyb9y4YRHNuHDG//m6ZHny5LGsV8bXLOP/OXPmlD2dUqwdPV7jX974mml8gwTbxHcNXbduHftQntf2suxjHp02a9YsS34uOh44cEB2WW8yQijzhhrKaJUAxizXntHjmKVHm117wLM7EMo844XcIAACICCbAHtzQQIBEAABEPCCANuxkYdRif5PmDDBixpRxEgE2K6cQrdu3ZyeHx07dhTYemo+dZWJY+n19uzZ06e65BRmmwykt+f4nB8wYICcKpAHBEBA4wT0OGbp0WYlnwbTp093OTavXLlSyaZQFwiAAAiYigCZqrfoLAiAAAgoSABCmYIwDVYV//LWtWtXpy8wI0eOFFJTU33qLYuGESIiItLrnjJlik/1ySkMoUwOJeQBAf0S0OOYpUeblX6GQChTmijqAwEQAIE0Aph6KTv2DhlBAARAwJ4Apl7a88BZGgH25Y3efvttyyL9ViYhISGWNcpYhJn1ktePv/zyi2WNN2sFW7ZssWyIYD1X4xFTL9WgijpBQBsE9Dhm6dFmNbyNqZdqUEWdIAACIMAIQDEEARAAARDwjgAiyrzjZuRSycnJQqdOndKjvfjbbNasWYVt27Yp1u0RI0bY1X/58mXF6nZVESLKXJHBdRDQNwE9jll6tFmtZwkiytQii3pBAATMTgARZZBLQQAEQMBLAogo8xKcQYuxL2/UuXNnWrFiRXoP+Q6f/Hny8ssvp1/z5eDx48dUunRpy+YJvJ4sWbJQdHS0L1XKKouIMlmYkAkEdEVAj2OWHm1W80mBiDI16aJuEAABMxOAUGZm76PvIAACPhGAUOYTPkMV5l/eWCQZ/fDDD3b9WrhwIRUqVMjumrsT9uud5TZbx4z4znyJiYkUExNDjx49YtvX36UdO3bY7XBZtWpVOnTokLsqFbkHoUwRjKgEBDRDQI9jlh5tVtvhEMrUJoz6QQAEzEoAQplZPY9+gwAI+EwAQpnPCA1RARe0oqKiaM2aNX7vD18LjYtYaqeZM2fSf//7X9FmevToQQsWLBC9h4sgAALaI6DHMUuPNvvD8xDK/EEZbYAACJiRAIQyM3odfQYBEFCEAIQyRTDquhIe8dWhQwdav359QPoxadIkGjJkiOptDx06lCZPnizaToMGDSyRbqI3cREEQEBTBPQ4ZunRZn85HUKZv0ijHRAAAbMRgFBmNo+jvyAAAooRgFCmGEpdVsS/vLVr1442bNgQMPvZhhLUvHlzVdvfuXMntWjRghISEly2wyPO+vbtS0FBQS7z4AYIgEBgCehxzNKjzf70MoQyf9JGWyAAAmYiAKHMTN5GX0EABBQlAKFMUZy6qoyvJcZFskBMt7QFdfHiRSpSpIjtJdnHt2/fppUrV6bn5+v/8OlNXBC7d+8eXbt2jY4fP05nz55Nz+PuIDIykipXrkyFCxem3LlzU0REBIWGhlJISEi6gFa/fn0qW7asu2pwDwRAQAUCehyz9GizCq5zWyWEMrd4cBMEQAAEvCYAocxrdCgIAiBgdgIQysz7DOCL6jdq1CigAJ555hniu2B6G8V15MgRqlKlil/7MH/+fOrZs6df20RjIAACZJkerbcxywjjrNrPPQhlahNG/SAAAmYlAKHMrJ5Hv0EABHwmAKHMZ4S6rWD58uX0xhtvBNT+SpUqERe7vE0Qyrwlh3IgoD8Cehyz9Gizv58ZEMr8TRztgQAImIUAhDKzeBr9BAEQUJwAhDLFkeqmwu3bt1Pjxo0Dam+XLl1o8eLFXttw/fp1qlChAj158oTCw8Mt/xkzZqSwsLD0fz51kv/za3wKJf/PkCGDpc2UlBSyTtfkUzb5WkLWR37M//k0Tut/XFwczZ07FxFlXnsMBUHAewJ6HLP0aLP3HvKuJIQy77ihFAiAAAhIEYBQJkUI90EABEDABQEIZS7A4DIIgAAIgAAIgIDqBCCUqY4YDYAACJiUAIQykzoe3QYBEPCdAIQy3xmiBhAAARAAARAAAe8IQCjzjhtKgQAIgIAUAQhlUoRwHwRAAARcEIBQ5gIMLoMACIAACIAACKhOAEKZ6ojRAAiAgEkJQCgzqePRbRAAAd8JQCjznSFqAAEQAAEQAAEQ8I4AhDLvuKEUCIAACEgRgFAmRQj3QQAEQMAFAQhlLsDgMgiAAAiAAAiAgOoEIJSpjhgNgAAImJQAhDKTOh7dBgEQ8J0AhDLfGaIGEAABEAABEAAB7whAKPOOG0qBAAiAgBQBCGVShHAfBEAABFwQgFDmAgwugwAIgAAIgAAIqE4AQpnqiNEACICASQlAKDOp49FtEAAB3wlAKPOdIWoAARAAARAAARDwjgCEMu+4oRQIgAAISBGAUCZFCPdBAARAwAUBCGUuwOAyCIAACIAACICA6gQglKmOGA2AAAiYlACEMpM6Ht0GARDwnQCEMt8ZogYQAAEQAAEQAAHvCEAo844bSoEACICAFAEIZVKEcB8EQAAEXBCAUOYCDC6DAAiAAAiAAAioTgBCmeqI0QAIgIBJCUAoM6nj0W0QAAHfCUAo850hagABEAABEAABEPCOAIQy77ihFAiAAAhIEYBQJkUI90EABEDABQEIZS7A4DIIgAAIgAAIgIDqBCCUqY4YDYAACJiUAIQykzoe3QYBEPCdAIQy3xmiBhAAARAAARAAAe8IQCjzjhtKgQAIgIAUAQhlUoRwHwRAAARcEIBQ5gIMLoMACIAACIAACKhOAEKZ6ojRAAiAgEkJQCgzqePRbRAAAd8JQCjznSFqAAEQAAEQAAEQ8I4AhDLvuKEUCIAACEgRgFAmRQj3QQAEQMAFAXdCWVRUFLVq1cpFSd8uV69enQoXLuxbJSgNAiAAAiAAAiCgOoHdu3fTzZs3VWln69at9O2334rWvXLlSmrfvr3oPVwEARAAARBwTwBCmXs+uAsCIAACLgm4E8pcFlLgRv369Wnnzp0K1IQqQAAEQAAEQAAE1CTQqFEj2rFjh5pNiNYNoUwUCy6CAAiAgCwCEMpkYUImEAABEHAmECihjEeUHThwwNkgXAEBEAABEAABENAUAQhlmnIHjAEBEAABWQQglMnChEwgAAIg4EwAQpkzE1wBARAAARAAARB4SgBC2VMWOAIBEAABvRCAUKYXT8FOEAABzRGAUKY5l8AgEAABEAABENAUAQhlmnIHjAEBEAABWQQglMnChEwgAAIg4EwAQpkzE1wBARAAARAAARB4SgBC2VMWOAIBEAABvRCAUKYXT8FOEAABzRGAUKY5l8AgEAABEAABENAUgebNm9PmzZv9bhMW8/c7cjQIAiBgIAIQygzkTHQFBEDAvwSio6OJfxBNTEz0a8P58+enNm3a+LVNNAYCIAACIAACIOA5gXPnztG2bds8L+hjiSZNmlDx4sV9rAXFQQAEQMCcBCCUmdPv6DUIgAAIgAAIgAAIgAAIgAAIgAAIgAAIgIADAQhlDkBwCgIgAAIgAAIgAAIgAAIgAAIgAAIgAAIgYE4CEMrM6Xf0GgRAAARAAARAAARAAARAAARAAARAAARAwIEAhDIHIDgFARAAARAAARAAARAAARAAARAAARAAARAwJwEIZeb0O3oNAiAAAiAAAiAAAiAAAiAAAiAAAiAAAiDgQABCmQMQnIIACIAACIAACIAACIAACIAACIAACIAACJiTAIQyc/odvQYBEAABEAABEAABEAABEAABEAABEAABEHAgAKHMAQhOQQAEQAAEQAAEQAAEQAAEQAAEQAAEQAAEzEkAQpk5/Y5egwAIgAAIgAAIgAAIgAAIgAAIgAAIgAAIOBCAUOYABKcgAAIgAAIgAAIgAAIgAAIgAAIgAAIgAALmJAChzJx+R69BAARAAARAAARAAARAAARAAARAAARAAAQcCEAocwCCUxAAARAAARAAARAAARAAARAAARAAARAAAXMSgFBmTr+j1yAAAiAAAiAAAiAAAiAAAiAAAiAAAiAAAg4EIJQ5AMEpCIAACIAACPxfO3ZIAAAAgDCsf2s8FT6JZDgIECBAgAABAgQIECDQFHCUNXfXmgABAgQIECBAgAABAgQIECBA4AQcZQciEiBAgAABAgQIECBAgAABAgQINAUcZc3dtSZAgAABAgQIECBAgAABAgQIEDgBR9mBiAQIECBAgAABAgQIECBAgAABAk0BR1lzd60JECBAgAABAgQIECBAgAABAgROwFF2ICIBAgQIECBAgAABAgQIECBAgEBTwFHW3F1rAgQIECBAgAABAgQIECBAgACBE3CUHYhIgAABAgQIECBAgAABAgQIECDQFHCUNXfXmgABAgQIECBAgAABAgQIECBA4AQcZQciEiBAgAABAgQIECBAgAABAgQINAUcZc3dtSZAgAABAgQIECBAgAABAgQIEDgBR9mBiAQIECBAgAABAgQIECBAgAABAk0BR1lzd60JECBAgAABAgQIECBAgAABAgROwFF2ICIBAgQIECBAgAABAgQIECBAgEBTwFHW3F1rAgQIECBAgAABAgQIECBAgACBE3CUHYhIgAABAgQIECBAgAABAgQIECDQFHCUNXfXmgABAgQIECBAgAABAgQIECBA4AQcZQciEiBAgAABAgQIECBAgAABAgQINAUcZc3dtSZAgAABAgQIECBAgAABAgQIEDgBR9mBiAQIECBAgAABAgQIECBAgAABAk0BR1lzd60JECBAgAABAgQIECBAgAABAgROwFF2ICIBAgQIECBAgAABAgQIECBAgEBTwFHW3F1rAgQIECBAgAABAgQIECBAgACBE3CUHYhIgAABAgQIECBAgAABAgQIECDQFHCUNXfXmgABAgQIECBAgAABAgQIECBA4AQcZQciEiBAgAABAgQIECBAgAABAgQINAUcZc3dtSZAgAABAgQIECBAgAABAgQIEDgBR9mBiAQIECBAgAABAgQIECBAgAABAk0BR1lzd60JECBAgAABAgQIECBAgAABAgROwFF2ICIBAgQIECBAgAABAgQIECBAgEBTwFHW3F1rAgQIECBAgAABAgQIECBAgACBE3CUHYhIgAABAgQIECBAgAABAgQIECDQFHCUNXfXmgABAgQIECBAgAABAgQIECBA4AQcZQciEiBAgAABAgQIECBAgAABAgQINAUcZc3dtSZAgAABAgQIECBAgAABAgQIEDgBR9mBiAQIECBAgAABAgQIECBAgAABAk0BR1lzd60JECBAgAABAgQIECBAgAABAgROwFF2ICIBAgQIECBAgAABAgQIECBAgEBTwFHW3F1rAgQIECBAgAABAgQIECBAgACBE3CUHYhIgAABAgQIECBAgAABAgQIECDQFHCUNXfXmgABAgQIECBAgAABAgQIECBA4AQcZQciEiBAgAABAgQIECBAgAABAgQINAUcZc3dtSZAgAABAgQIECBAgAABAgQIEDgBR9mBiAQIECBAgAABAgQIECBAgAABAk0BR1lzd60JECBAgAABAgQIECBAgAABAgROYAmBS0HCGRr1AAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "0b3375d5",
   "metadata": {},
   "source": [
    "$$\\require{nicematrix}$$\n",
    "\n",
    "# 2. The LP steady-state optimizer in DMC\n",
    "\n",
    "## MVs and CVs\n",
    "\n",
    "In multivariable controllers, we have Manipulated Variables (MVs) which are handles that can be adjusted to move the process. The effects of moving the MVs are reflected in the Controlled Variable (CVs).\n",
    "\n",
    "For example, adjusting a reboiler's steam flow rate (MV) will affect the temperature of the column bottoms (CV).\n",
    "\n",
    "## Gain Matrix\n",
    "\n",
    "The steady-state relationship between an MV and a CV is captured by its **steady-state process gain**. In multivariable controllers, we have a **gain matrix** that captures the relationships between multiple MVs and CVs.\n",
    "\n",
    "Consider a $2\\times2$ gain matrix, $G$ with 2 MVs and 2 CVs. There will be 4 elements, where each $G_{ij}$ describes the relationship between the $i$-th CV and $j$-th MV.\n",
    "\n",
    "$$\n",
    "G =\n",
    "    \\begin{bmatrix}\n",
    "        G_{11} & G_{12} \\\\\n",
    "        G_{21} & G_{22}\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "> __Note__ that mathematically, the CVs are rows, and MVs are columns. This is the opposite of the usual convention in commercial MPC software, where the MVs are plotted as rows, and CVs as columns.\n",
    "\n",
    "<div>\n",
    "<br>\n",
    "<img src=\"attachment:image-3.png\" width=\"300\"/>\n",
    "</div>\n",
    "\n",
    "The equation relating the CVs and MVs through the gain matrix is given by:\n",
    "\n",
    "$$ \\Delta CV = G \\cdot \\Delta MV $$\n",
    "\n",
    "Multiplying the terms in the matrix, we get:\n",
    "\n",
    "$$\n",
    "\\Delta \\text{CV}_{1} = G_{11} \\Delta \\text{MV}_{1} + G_{12} \\Delta \\text{MV}_{2} \\\\ \n",
    "\\Delta \\text{CV}_{2} = G_{21} \\Delta \\text{MV}_{1} + G_{22} \\Delta \\text{MV}_{2} \n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2c321ce",
   "metadata": {},
   "source": [
    "\n",
    "##  2.1 The debutanizer process\n",
    "\n",
    "Let's look at a debutanizer process, which separates C4- molecules at the top of the column from heavier components at the bottom.\n",
    "\n",
    "<div>\n",
    "<br>\n",
    "<img src=\"./assets/img/debut.png\" width=\"600\"/>\n",
    "</div>\n",
    "\n",
    "From: *O'Connor, D. L., Siang, L. C., & Elnawawi, S. (2022). A Simple Discretization Scheme for Gain Matrix Conditioning. arXiv preprint arXiv:2202.13291.*\n",
    "\n",
    "To simplify the problem, let's just consider 2 MVs and 2 CVs.\n",
    "\n",
    "- **MV1:** Reboiler Temperature Control (TC-REBOILER)\n",
    "- **MV2:** Reflux Flow Rate Control (FC-REFLUX)\n",
    "- **CV1:** Bottoms RVP (AI-BTM-RVP)\n",
    "- **CV2:** Overhead C5 Content (AI-TOP-C5)\n",
    "\n",
    "The gain matrix is given as:\n",
    "$$\n",
    "G =\n",
    "    \\begin{bmatrix}\n",
    "        -0.200 & -0.072 \\\\\n",
    "        0.125 & -0.954\n",
    "    \\end{bmatrix}\n",
    "$$\n",
    "\n",
    "The gain matrix tells us that:\n",
    "- Every unit increase in reboiler temperature reduces RVP by 0.20 units and increases C5 by 0.125 units.\n",
    "- This makes sense, if we fire the reboiler harder, we boil stuff up the top so the bottom RVP decreases (less volatile, more heavier components), and some of the heavier components go up the column, so overhead C5 increases.\n",
    "- Every unit increase in reflux reduces RVP by 0.072 units and reduces C5 by 0.954 units.\n",
    "- This also makes sense, because increasing reflux improves separation (up to a certain point, considering flooding limits etc.)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "51f42715",
   "metadata": {},
   "source": [
    "## 2.2 MV-CV equations\n",
    "\n",
    "Using the gain matrix, the CV relationship can be written in terms of its MVs, starting with:\n",
    "\n",
    "$$\n",
    "\\Delta \\text{CV}_{1} = G_{11} \\Delta \\text{MV}_{1} + G_{12} \\Delta \\text{MV}_{2} \\\\ \n",
    "\\Delta \\text{CV}_{2} = G_{21} \\Delta \\text{MV}_{1} + G_{22} \\Delta \\text{MV}_{2} \n",
    "$$\n",
    "\n",
    "We can impose upper and lower limits on the MVs:\n",
    "\n",
    "$$\n",
    "\\text{MV}_{1, \\text{Lo}} \\leq \\text{MV}_{1} \\leq \\text{MV}_{1, \\text{Hi}}\\\\ \n",
    "\\text{MV}_{2, \\text{Lo}} \\leq \\text{MV}_{2} \\leq \\text{MV}_{2, \\text{Hi}}\\\\\n",
    "$$\n",
    "\n",
    "As well as the CVs:\n",
    "\n",
    "$$\n",
    "\\text{CV}_{1, \\text{Lo}} \\leq \\text{CV}_{1} \\leq \\text{CV}_{1, \\text{Hi}}\\\\ \n",
    "\\text{CV}_{2, \\text{Lo}} \\leq \\text{CV}_{2} \\leq \\text{CV}_{2, \\text{Hi}}\\\\\n",
    "$$\n",
    "\n",
    "Since the CVs are related to the MVs by the gain matrix, we can substitute the equations to get CV limits in terms of MV movements:\n",
    "\n",
    "$$\n",
    "G_{11} \\Delta \\text{MV}_{1} + G_{12} \\Delta \\text{MV}_{2} \\leq \\Delta \\text{CV}_{1, \\text{Hi}}\\\\\n",
    "G_{11} \\Delta \\text{MV}_{1} + G_{12} \\Delta \\text{MV}_{2} \\geq \\Delta \\text{CV}_{1, \\text{Lo}}\\\\\n",
    "G_{21} \\Delta \\text{MV}_{1} + G_{22} \\Delta \\text{MV}_{2} \\leq \\Delta \\text{CV}_{2, \\text{Hi}} \\\\\n",
    "G_{21} \\Delta \\text{MV}_{1} + G_{22} \\Delta \\text{MV}_{2} \\geq \\Delta \\text{CV}_{2, \\text{Lo}} \\\\\n",
    "$$\n",
    "\n",
    "Let's see what this looks like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5d2d53d8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Just importing libraries and tweaking the plot settings\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib widget\n",
    "from matplotlib.ticker import StrMethodFormatter\n",
    "from IPython.display import display, clear_output\n",
    "\n",
    "fsize = 15\n",
    "tsize = 10\n",
    "tdir = 'in'\n",
    "major = 5.0\n",
    "minor = 3.0\n",
    "lwidth = 0.8\n",
    "lhandle = 2.0\n",
    "plt.style.use('default')\n",
    "plt.rcParams['font.size'] = fsize\n",
    "plt.rcParams['legend.fontsize'] = tsize\n",
    "plt.rcParams['xtick.direction'] = tdir\n",
    "plt.rcParams['ytick.direction'] = tdir\n",
    "plt.rcParams['xtick.major.size'] = major\n",
    "plt.rcParams['xtick.minor.size'] = minor\n",
    "plt.rcParams['ytick.major.size'] = 5.0\n",
    "plt.rcParams['ytick.minor.size'] = 3.0\n",
    "plt.rcParams['axes.linewidth'] = lwidth\n",
    "plt.rcParams['legend.handlelength'] = lhandle"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30414a69",
   "metadata": {},
   "source": [
    "# Gains and CV Limits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58065fad",
   "metadata": {},
   "outputs": [],
   "source": [
    "G11 = -0.200\n",
    "G12 = -0.072\n",
    "G21 = 0.125\n",
    "G22 = -0.954\n",
    "\n",
    "CV1Lo = -1\n",
    "CV1Hi = 1\n",
    "\n",
    "CV2Lo = -2\n",
    "CV2Hi = 2.5\n",
    "\n",
    "# standard figure size\n",
    "fig_std = (8, 5)\n",
    "\n",
    "# standard figure parameters (constraints)\n",
    "limits = 10\n",
    "d = np.linspace(-limits, limits, 100)\n",
    "x,y = np.meshgrid(d,d)\n",
    "\n",
    "# constraint formulation in terms of x and y\n",
    "c1 = G11*x+G12*y <= CV1Hi\n",
    "c2 = G11*x+G12*y >= CV1Lo\n",
    "c3 = G21*x+G22*y <= CV2Hi\n",
    "c4 = G21*x+G22*y >= CV2Lo\n",
    "\n",
    "# equation of a line, y = mx + c\n",
    "y_c1 = (CV1Hi - G11*d)/G12\n",
    "y_c2 = (CV1Lo - G11*d)/G12\n",
    "y_c3 = (CV2Hi - G21*d)/G22\n",
    "y_c4 = (CV2Lo - G21*d)/G22"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc30de46",
   "metadata": {},
   "source": [
    "## 2.3 The Feasible Region\n",
    "\n",
    "The shaded area below is known as the `feasible region` or `solution space` in the LP problem. It is defined based on the current process conditions and the distance of each variable from its constraints."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70b9e060",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=fig_std)\n",
    "plt.imshow((c1 & c2 & c3 & c4).astype(int), extent=(x.min(),x.max(),y.min(),y.max()), origin=\"lower\", cmap=\"Blues\", alpha=0.1);\n",
    "plt.plot(d, y_c1, '-r', label='CV1_Hi');\n",
    "plt.plot(d, y_c2, '--r', label='CV1_Lo');\n",
    "plt.plot(d, y_c3, '-b', label='CV2_Hi');\n",
    "plt.plot(d, y_c4, '--b', label='CV2_Lo');\n",
    "plt.plot(0,0,'kx');\n",
    "\n",
    "plt.xlim((-limits, limits))\n",
    "plt.ylim((-limits, limits))\n",
    "plt.xlabel(r'$\\Delta MV_1$')\n",
    "plt.ylabel(r'$\\Delta MV_2$')\n",
    "plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}')); # No decimal places\n",
    "plt.legend(loc='lower center', bbox_to_anchor=(0.5, 1.03), ncol=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b4adbbfa",
   "metadata": {},
   "source": [
    "## 2.4 The optimum point (LP solution)\n",
    "The feasible region tells us that the LP optimizer is allowed to move $\\Delta MV_1$ and $\\Delta MV_2$ within the shaded region to honour the CV limits.\n",
    "\n",
    "The question now is, out of all the possible points in the feasible region, which one should the optimizer pick and why? Recall in the previous chapter that for the cat food production problem, the objective function is defined as finding combinations of chicken, $x_1$ and beef, $x_2$, given the cost of chicken $c_1$, and cost of beef $c_2$, such that the total production costs are minimized:\n",
    "\n",
    "$$\\min_{x_1, x_2} f(x_1, x_2) = c_1 x_1 + c_2 x_2$$\n",
    "\n",
    "The objective function in DMC is defined in a similar manner. The equations below are a simplified version of the actual objective function (see *Sorensen, R. C., & Cutler, C. R. (1998)* for details).\n",
    "\n",
    "We want to assign an 'LP cost' to each MV, based on the economics and desired directionality of MV movement. For 2 MVs, we get:\n",
    "\n",
    "$$\\min_{\\Delta MV_1, \\Delta MV_2} f(\\Delta MV_1, \\Delta MV_2) = c_1 \\Delta MV_1 + c_2 \\Delta MV_2$$\n",
    "\n",
    "As a rule of thumb, a negative LP cost would incentivize the DMC LP to maximize that variable, and likewise, a positive cost would incentivize the DMC LP to minimize that variable. However, there are exceptions as we will see later on.\n",
    "\n",
    "For now, let's assume that we have the following LP costs:\n",
    "\n",
    "- $c_1 = -1$\n",
    "- $c_2 = -1$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6915ed5e",
   "metadata": {},
   "outputs": [],
   "source": [
    "cost_MV1 = -1\n",
    "cost_MV2 = -1"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "775599bb",
   "metadata": {},
   "source": [
    "## 3 Set up the LP problem in PuLP\n",
    "\n",
    "Combining the previous sections, our DMC LP problem is defined as follows.\n",
    "\n",
    "### 3.1. Decision Variables\n",
    "The **decision variables** are our MVs.\n",
    "\n",
    "- **MV1:** Reboiler Temperature Control (TC-REBOILER)\n",
    "- **MV2:** Reflux Flow Rate Control (FC-REFLUX)\n",
    "\n",
    "### 3.2. Objective Function\n",
    "The **objective function** is defined as minimizing the cost of moving our MVs:\n",
    "\n",
    "$$\\min_{\\Delta MV_1, \\Delta MV_2} f(\\Delta MV_1, \\Delta MV_2) = c_1 \\Delta MV_1 + c_2 \\Delta MV_2$$\n",
    "\n",
    "### 3.3 Constraints\n",
    "Our 2x2 DMC problem is subjected to these constraints\n",
    "\n",
    "#### 3.3.1 MV Constraints\n",
    "\n",
    "$$\n",
    "\\Delta \\text{MV}_{1} \\leq \\Delta \\text{MV}_{1, \\text{Hi}}\\\\ \n",
    "\\Delta \\text{MV}_{1} \\geq \\Delta \\text{MV}_{1, \\text{Lo}}\\\\\n",
    "\\Delta \\text{MV}_{2} \\leq \\Delta \\text{MV}_{2, \\text{Hi}}\\\\\n",
    "\\Delta \\text{MV}_{2} \\geq \\Delta \\text{MV}_{2, \\text{Lo}} \\\\\n",
    "$$\n",
    "\n",
    "#### 3.3.2 CV Constraints (written in terms of MVs)\n",
    "\n",
    "- **CV1:** Bottoms RVP (AI-BTM-RVP)\n",
    "\n",
    "$$\n",
    "G_{11} \\Delta \\text{MV}_{1} + G_{12} \\Delta \\text{MV}_{2} \\leq \\Delta \\text{CV}_{1, \\text{Hi}}\\\\\n",
    "G_{11} \\Delta \\text{MV}_{1} + G_{12} \\Delta \\text{MV}_{2} \\geq \\Delta \\text{CV}_{1, \\text{Lo}}\\\\\n",
    "$$\n",
    "\n",
    "- **CV2:** Overhead C5 Content (AI-TOP-C5)\n",
    "\n",
    "$$\n",
    "G_{21} \\Delta \\text{MV}_{1} + G_{22} \\Delta \\text{MV}_{2} \\leq \\Delta \\text{CV}_{2, \\text{Hi}} \\\\\n",
    "G_{21} \\Delta \\text{MV}_{1} + G_{22} \\Delta \\text{MV}_{2} \\geq \\Delta \\text{CV}_{2, \\text{Lo}} \\\\\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "656132e3",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import PuLP modeler functions\n",
    "from pulp import *"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e949f7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The 2 variables MV1 and MV2 are created with a lower limit of -10\n",
    "MV1=LpVariable(\"MV1\",-limits)\n",
    "MV2=LpVariable(\"MV2\",-limits)\n",
    "\n",
    "# Create the 'prob' variable to contain the problem data\n",
    "prob = LpProblem(\"DMC_problem\",LpMinimize)\n",
    "\n",
    "# The objective function is added to 'prob' first\n",
    "prob += cost_MV1*MV1 + cost_MV2*MV2, \"Cost function of MVs\"\n",
    "\n",
    "# constraint formulation in terms of MV1 and MV2\n",
    "prob += G11*MV1+G12*MV2 <= CV1Hi, \"CV1 High Limit\"\n",
    "prob += G11*MV1+G12*MV2 >= CV1Lo, \"CV1 Low Limit\"\n",
    "prob += G21*MV1+G22*MV2 <= CV2Hi, \"CV2 High Limit\"\n",
    "prob += G21*MV1+G22*MV2 >= CV2Lo, \"CV2 Low Limit\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1e9d6586",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39d32d81",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The problem is solved using PuLP's choice of Solver\n",
    "prob.solve()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e7768ac6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Each of the variables is printed with it's resolved optimum value\n",
    "for v in prob.variables():\n",
    "    print(v.name, \"=\", v.varValue)\n",
    "    \n",
    "soln = [v.varValue for v in prob.variables()]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5263ab1a",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=fig_std)\n",
    "plt.imshow((c1 & c2 & c3 & c4).astype(int), extent=(x.min(),x.max(),y.min(),y.max()), origin=\"lower\", cmap=\"Blues\", alpha=0.1);\n",
    "plt.plot(d, y_c1, '-r', label='CV1_Hi');\n",
    "plt.plot(d, y_c2, '--r', label='CV1_Lo');\n",
    "plt.plot(d, y_c3, '-b', label='CV2_Hi');\n",
    "plt.plot(d, y_c4, '--b', label='CV2_Lo');\n",
    "plt.plot(0,0,'kx');\n",
    "plt.plot(soln[0], soln[1], 'ok');\n",
    "\n",
    "plt.xlim((-limits, limits))\n",
    "plt.ylim((-limits, limits))\n",
    "plt.xlabel(r'$\\Delta MV_1$')\n",
    "plt.ylabel(r'$\\Delta MV_2$')\n",
    "plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}')); # No decimal places\n",
    "plt.legend(loc='lower center', bbox_to_anchor=(0.5, 1.03), ncol=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d803a51e",
   "metadata": {},
   "source": [
    "---------------------------\n",
    "## Objective Function"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "58a6b08a",
   "metadata": {},
   "source": [
    "The objective function value, $V$, can be accessed using"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "487ad44e",
   "metadata": {},
   "outputs": [],
   "source": [
    "value(prob.objective)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d46b9cc",
   "metadata": {},
   "source": [
    "Recall that the objective function is:\n",
    "\n",
    "$$\\min_{\\Delta MV_1, \\Delta MV_2} f(\\Delta MV_1, \\Delta MV_2) = c_1 \\Delta MV_1 + c_2 \\Delta MV_2$$\n",
    "\n",
    "So our line would be:\n",
    "\n",
    "$$c_1 x + c_2 y = V$$\n",
    "\n",
    "Rearrange the equation to get:\n",
    "\n",
    "$$ y = \\frac{1}{c_2} V - \\frac{c_1}{c_2} x $$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2dc475e",
   "metadata": {},
   "source": [
    "Let's plot the line defining the objective function at the optimum solution. This plot will also include a vector field showing the direction of optimization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b3e6805",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# grid space for vector definition\n",
    "dvec = np.linspace(-limits + (0.1*limits), limits - (0.1*limits), 12)\n",
    "xv, yv = np.meshgrid(dvec, dvec)\n",
    "\n",
    "V = value(prob.objective)\n",
    "y_obj = (1/cost_MV2)*V - (cost_MV1/cost_MV2)*d\n",
    "z_obj = (cost_MV1 * xv) + (cost_MV2 * yv)\n",
    "\n",
    "plt.figure(figsize=fig_std)\n",
    "plt.imshow((c1 & c2 & c3 & c4).astype(int), extent=(x.min(),x.max(),y.min(),y.max()), origin=\"lower\", cmap=\"Blues\", alpha=0.1);\n",
    "plt.plot(d, y_c1, '-r', label='CV1_Hi');\n",
    "plt.plot(d, y_c2, '--r', label='CV1_Lo');\n",
    "plt.plot(d, y_c3, '-b', label='CV2_Hi');\n",
    "plt.plot(d, y_c4, '--b', label='CV2_Lo');\n",
    "plt.plot(0,0,'kx');\n",
    "plt.plot(soln[0], soln[1], 'ok');\n",
    "\n",
    "# the obj func\n",
    "plt.plot(d, y_obj, '-.k');\n",
    "\n",
    "# plot vector field to show the direction of optimization (direction of decreasing cost)\n",
    "plt.quiver(xv, yv,-cost_MV1,-cost_MV2, z_obj, cmap='gray', headwidth=4, width=0.003, scale=40, alpha=0.5)\n",
    "\n",
    "# plot parameters\n",
    "plt.xlim((-limits, limits))\n",
    "plt.ylim((-limits, limits))\n",
    "plt.xlabel(r'$\\Delta MV_1$')\n",
    "plt.ylabel(r'$\\Delta MV_2$')\n",
    "plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}')); # No decimal places\n",
    "plt.legend(loc='lower center', bbox_to_anchor=(0.5, 1.03), ncol=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8c8866b8",
   "metadata": {},
   "source": [
    "#### A note on the vector field\n",
    "> Note: the vector field in the preceding plot shows the direction of *decreasing* cost, so this is the direction of the feasible region that the LP will try to drive the solution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55ab981c",
   "metadata": {},
   "source": [
    "### Create a function to run the LP solver easily\n",
    "\n",
    "This will define a ```run_pulp()``` function, which contains the previously used LP solver. This way, we can simply call on ```run_pulp()``` instead of writing out all the code we used to obtain the LP solution. Writing it as a function will also make it easier to modify costs and visualize how the LP solution will change."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "110b9543",
   "metadata": {},
   "outputs": [],
   "source": [
    "def run_pulp(cost_MV1, cost_MV2, MV1_Lo=-limits, MV2_Lo=-limits, MV1_Hi=limits, MV2_Hi=limits):\n",
    "    '''\n",
    "    Run pulp solver for LP problem with the same gain matrix and CV limits as the example. Returns a tuple with:\n",
    "    (LP solution array [MV1_solution, MV2_solution], objective_value)\n",
    "    '''\n",
    "    \n",
    "    # Gains and CV limits\n",
    "    G11 = -0.200\n",
    "    G12 = -0.072\n",
    "    G21 = 0.125\n",
    "    G22 = -0.954\n",
    "\n",
    "    CV1Lo = -1\n",
    "    CV1Hi = 1\n",
    "\n",
    "    CV2Lo = -2\n",
    "    CV2Hi = 2.5\n",
    "\n",
    "    prob = LpProblem(\"DMC_2\", LpMinimize)\n",
    "    \n",
    "    MV1 = LpVariable(\"MV1\", MV1_Lo)\n",
    "    MV2 = LpVariable(\"MV2\", MV2_Lo)\n",
    "    \n",
    "    prob += cost_MV1 * MV1 + cost_MV2 * MV2\n",
    "    \n",
    "    # constraint formulation in terms of MV1 and MV2\n",
    "    prob += G11*MV1+G12*MV2 <= CV1Hi, \"CV1 High Limit\"\n",
    "    prob += G11*MV1+G12*MV2 >= CV1Lo, \"CV1 Low Limit\"\n",
    "    prob += G21*MV1+G22*MV2 <= CV2Hi, \"CV2 High Limit\"\n",
    "    prob += G21*MV1+G22*MV2 >= CV2Lo, \"CV2 Low Limit\"\n",
    "    prob += MV1 >= MV1_Lo, \"MV1 Low Limit\"\n",
    "    prob += MV1 <= MV1_Hi, \"MV1 High Limit\"\n",
    "    prob += MV2 >= MV2_Lo, \"MV2 Low Limit\"\n",
    "    prob += MV2 <= MV2_Hi, \"MV2 High Limit\"\n",
    "    \n",
    "    # The problem is solved using PuLP's choice of Solver\n",
    "    prob.solve(PULP_CBC_CMD(msg=0))\n",
    "    \n",
    "    # Each of the variables is printed with it's resolved optimum value\n",
    "#     for v in prob.variables():\n",
    "#         print(v.name, \"=\", v.varValue)\n",
    "\n",
    "    soln = [v.varValue for v in prob.variables()]\n",
    "    \n",
    "    obj_value = value(prob.objective)\n",
    "    \n",
    "    return (soln, obj_value)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "07124356",
   "metadata": {},
   "source": [
    "## Changes in the objective function\n",
    "\n",
    "As mentioned previously, a rule of thumb is that variables with negative LP costs are maximized and vice versa. However, the direction DMC will optimize a variable will be determined by the solution to the objective function.\n",
    "\n",
    "Recall that the objective function is:\n",
    "\n",
    "$$\\min_{\\Delta MV_1, \\Delta MV_2} f(\\Delta MV_1, \\Delta MV_2) = c_1 \\Delta MV_1 + c_2 \\Delta MV_2$$\n",
    "\n",
    "### One negative, one positive cost\n",
    "\n",
    "Let's see how changing the LP costs will affect the solution. Assume that we changed the costs as follows:\n",
    "\n",
    "- $c_1 = -1$\n",
    "- $c_2 = +1$ \\(same as before)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b842ad4e",
   "metadata": {},
   "outputs": [],
   "source": [
    "lpc_MV1 = -1\n",
    "lpc_MV2 = 1\n",
    "\n",
    "soln, V = run_pulp(lpc_MV1, lpc_MV2)\n",
    "y_obj = (1/lpc_MV2)*V - (lpc_MV1/lpc_MV2)*d\n",
    "z_obj = (lpc_MV1 * xv) + (lpc_MV2 * yv)\n",
    "\n",
    "plt.figure(figsize=fig_std)\n",
    "plt.imshow((c1 & c2 & c3 & c4).astype(int), extent=(x.min(),x.max(),y.min(),y.max()), origin=\"lower\", cmap=\"Blues\", alpha=0.1);\n",
    "plt.plot(d, y_c1, '-r', label='CV1_Hi');\n",
    "plt.plot(d, y_c2, '--r', label='CV1_Lo');\n",
    "plt.plot(d, y_c3, '-b', label='CV2_Hi');\n",
    "plt.plot(d, y_c4, '--b', label='CV2_Lo');\n",
    "\n",
    "plt.plot(0,0,'kx');\n",
    "# the soln\n",
    "plt.plot(soln[0], soln[1], 'ok');\n",
    "\n",
    "# the obj func\n",
    "plt.plot(d, y_obj, '-.k');\n",
    "\n",
    "# plot vector field to show the direction of optimization (direction of decreasing cost)\n",
    "plt.quiver(xv, yv,-lpc_MV1,-lpc_MV2, z_obj, cmap='gray', headwidth=4, width=0.003, scale=40, alpha=0.5)\n",
    "\n",
    "plt.xlim((-limits, limits))\n",
    "plt.ylim((-limits, limits))\n",
    "plt.xlabel(r'$\\Delta MV_1$')\n",
    "plt.ylabel(r'$\\Delta MV_2$')\n",
    "plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}')); # No decimal places\n",
    "plt.legend(loc='lower center', bbox_to_anchor=(0.5, 1.03), ncol=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "33da56cb",
   "metadata": {},
   "source": [
    "Notice that MV2 is now being minimized, while MV1 is still being maximized. This is because the LP cost is positive, meaning that increasing MV2 will incur additional costs for the objective function.\n",
    "\n",
    "### Both positive costs\n",
    "\n",
    "Similarly, if we set $c_1 = +1$ and $c_2 = +1$, we should see that both variables will now be minimized."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2876471",
   "metadata": {},
   "outputs": [],
   "source": [
    "lpc_MV1 = 1\n",
    "lpc_MV2 = 1\n",
    "\n",
    "soln, V = run_pulp(lpc_MV1, lpc_MV2)\n",
    "y_obj = (1/lpc_MV2)*V - (lpc_MV1/lpc_MV2)*d\n",
    "z_obj = (lpc_MV1 * xv) + (lpc_MV2 * yv)\n",
    "\n",
    "plt.figure(figsize=fig_std)\n",
    "plt.imshow((c1 & c2 & c3 & c4).astype(int), extent=(x.min(),x.max(),y.min(),y.max()), origin=\"lower\", cmap=\"Blues\", alpha=0.1);\n",
    "plt.plot(d, y_c1, '-r', label='CV1_Hi');\n",
    "plt.plot(d, y_c2, '--r', label='CV1_Lo');\n",
    "plt.plot(d, y_c3, '-b', label='CV2_Hi');\n",
    "plt.plot(d, y_c4, '--b', label='CV2_Lo');\n",
    "\n",
    "plt.plot(0,0,'kx');\n",
    "plt.plot(soln[0], soln[1], 'ok');\n",
    "\n",
    "# the obj func\n",
    "plt.plot(d, y_obj, '-.k');\n",
    "\n",
    "# plot vector field to show the direction of optimization (direction of decreasing cost)\n",
    "plt.quiver(xv, yv,-lpc_MV1,-lpc_MV2, z_obj, cmap='gray', headwidth=4, width=0.003, scale=40, alpha=0.5)\n",
    "\n",
    "plt.xlim((-limits, limits))\n",
    "plt.ylim((-limits, limits))\n",
    "plt.xlabel(r'$\\Delta MV_1$')\n",
    "plt.ylabel(r'$\\Delta MV_2$')\n",
    "plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}')); # No decimal places\n",
    "plt.legend(loc='lower center', bbox_to_anchor=(0.5, 1.03), ncol=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "daac1e66",
   "metadata": {},
   "source": [
    "### Caution: costs do not guarantee the directions\n",
    "\n",
    "The rule of thumb that negative cost &rarr; maximized MV and positive cost &rarr; minimized MV does not always apply. That is only the 'preferred' direction, and not guaranteed.\n",
    "\n",
    "The direction that each variable will be driven is dependent on exactly where the LP solution lies, so some variables with negative LP cost may end up being minimized and vice versa.\n",
    "\n",
    "For example, if we set $c_2 = +0.3$ (i.e. where it is still positive):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "61e1f71d",
   "metadata": {},
   "outputs": [],
   "source": [
    "lpc_MV1 = 1\n",
    "lpc_MV2 = 0.3\n",
    "\n",
    "soln, V = run_pulp(lpc_MV1, lpc_MV2)\n",
    "y_obj = (1/lpc_MV2)*V - (lpc_MV1/lpc_MV2)*d\n",
    "z_obj = (lpc_MV1 * xv) + (lpc_MV2 * yv)\n",
    "\n",
    "plt.figure(figsize=fig_std)\n",
    "plt.imshow((c1 & c2 & c3 & c4).astype(int), extent=(x.min(),x.max(),y.min(),y.max()), origin=\"lower\", cmap=\"Blues\", alpha=0.1);\n",
    "plt.plot(d, y_c1, '-r', label='CV1_Hi');\n",
    "plt.plot(d, y_c2, '--r', label='CV1_Lo');\n",
    "plt.plot(d, y_c3, '-b', label='CV2_Hi');\n",
    "plt.plot(d, y_c4, '--b', label='CV2_Lo');\n",
    "\n",
    "plt.plot(0,0,'kx');\n",
    "plt.plot(soln[0], soln[1], 'ok');\n",
    "\n",
    "# the obj func\n",
    "plt.plot(d, y_obj, '-.k');\n",
    "\n",
    "# plot vector field to show the direction of optimization (direction of decreasing cost)\n",
    "plt.quiver(xv, yv,-lpc_MV1,-lpc_MV2, z_obj, cmap='gray', headwidth=4, width=0.003, scale=40, alpha=0.5)\n",
    "\n",
    "plt.xlim((-limits, limits))\n",
    "plt.ylim((-limits, limits))\n",
    "plt.xlabel(r'$\\Delta MV_1$')\n",
    "plt.ylabel(r'$\\Delta MV_2$')\n",
    "plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}')); # No decimal places\n",
    "plt.legend(loc='lower center', bbox_to_anchor=(0.5, 1.03), ncol=4);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2af23fae",
   "metadata": {},
   "source": [
    "#### Observe!\n",
    "\n",
    "> Even though both costs are positive, MV2 is being maximized. \n",
    "\n",
    "With both costs being positive, increasing either MV will increase the overall cost in the objective function. However, the costs incurred by increasing MV2 are not great enough to offset the profit generated by minimizing MV1."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c7b984f8",
   "metadata": {},
   "source": [
    "## Changes in MV limits\n",
    "\n",
    "All of these examples are assuming that we have constant CV and MV limits. In practice, we use the limits to move the variables around as necessary, so taking variable MV limits into account is an important aspect of the LP solution.\n",
    "\n",
    "Currently, both variables have limits between $(-10, 10)$. Let us see how the solution changes if we clamp the lower limit of MV1 at $-2.0$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2670205b",
   "metadata": {},
   "outputs": [],
   "source": [
    "lpc_MV1 = 1\n",
    "lpc_MV2 = 0.2\n",
    "\n",
    "MV1_lims = (-2.0, 10.0)\n",
    "MV2_lims = (-10.0, 10.0)\n",
    "\n",
    "m1 = (x >= MV1_lims[0]) & (x <= MV1_lims[1])\n",
    "m2 = (y >= MV2_lims[0]) & (y <= MV2_lims[1])\n",
    "\n",
    "soln, V = run_pulp(lpc_MV1, lpc_MV2, MV1_Lo=MV1_lims[0])\n",
    "y_obj = (1/lpc_MV2)*V - (lpc_MV1/lpc_MV2)*d\n",
    "z_obj = (lpc_MV1 * xv) + (lpc_MV2 * yv)\n",
    "\n",
    "plt.figure(figsize=fig_std)\n",
    "plt.imshow((c1 & c2 & c3 & c4 & m1 & m2).astype(int), extent=(x.min(),x.max(),y.min(),y.max()), origin=\"lower\", cmap=\"Blues\", alpha=0.1);\n",
    "plt.plot(d, y_c1, '-r', label='CV1_Hi');\n",
    "plt.plot(d, y_c2, '--r', label='CV1_Lo');\n",
    "plt.plot(d, y_c3, '-b', label='CV2_Hi');\n",
    "plt.plot(d, y_c4, '--b', label='CV2_Lo');\n",
    "\n",
    "plt.plot(MV1_lims[0] * np.ones(d.shape), d, 'g--', label='MV1_Lo');\n",
    "# plt.plot(MV1_lims[1] * np.ones(d.shape), d, 'g-', label='MV1_Hi');\n",
    "# plt.plot(MV2_lims[0] * np.ones(d.shape), d, 'g--', label='MV2_Lo');\n",
    "# plt.plot(MV2_lims[1] * np.ones(d.shape), d, 'g-', label='MV2_Hi');\n",
    "\n",
    "plt.plot(0,0,'kx');\n",
    "plt.plot(soln[0], soln[1], 'ok');\n",
    "\n",
    "# the obj func\n",
    "plt.plot(d, y_obj, '-.k');\n",
    "\n",
    "# plot vector field to show the direction of optimization (direction of decreasing cost)\n",
    "plt.quiver(xv, yv,-lpc_MV1,-lpc_MV2, z_obj, cmap='gray', headwidth=4, width=0.003, scale=40, alpha=0.5)\n",
    "\n",
    "plt.xlim((-limits, limits))\n",
    "plt.ylim((-limits, limits))\n",
    "plt.xlabel(r'$\\Delta MV_1$')\n",
    "plt.ylabel(r'$\\Delta MV_2$')\n",
    "plt.gca().yaxis.set_major_formatter(StrMethodFormatter('{x:,.0f}')); # No decimal places\n",
    "plt.legend(loc='lower center', bbox_to_anchor=(0.5, 1.02), ncol=5);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45cf5721",
   "metadata": {},
   "source": [
    "---------------------------------\n",
    "\n",
    "## Interactive LP solver\n",
    "\n",
    "Now let's look at an interactive example, where the LP costs and CV limits can be moved around to see how the LP solution is formulated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c127fa0c",
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "from ipywidgets.widgets.interaction import interact\n",
    "import ipywidgets.widgets as widgets\n",
    "from ipywidgets import Layout\n",
    "\n",
    "limits = 10\n",
    "d = np.linspace(-limits, limits, 100)\n",
    "x,y = np.meshgrid(d,d)\n",
    "stepsize = 0.1\n",
    "\n",
    "v_new = 0\n",
    "\n",
    "# mesh for vector field\n",
    "dvec = np.linspace(-limits + (0.1*limits), limits - (0.1*limits), 12)\n",
    "xv, yv = np.meshgrid(dvec, dvec)\n",
    "\n",
    "# constraint formulation in terms of x and y\n",
    "c1 = G11*x+G12*y <= CV1Hi\n",
    "c2 = G11*x+G12*y >= CV1Lo\n",
    "c3 = G21*x+G22*y <= CV2Hi\n",
    "c4 = G21*x+G22*y >= CV2Lo\n",
    "\n",
    "MV1_lims = (-9.0, 8)\n",
    "MV2_lims = (-9.0, 8)\n",
    "\n",
    "m1 = (x >= MV1_lims[0]) & (x <= MV1_lims[1])\n",
    "m2 = (y >= MV2_lims[0]) & (y <= MV2_lims[1])\n",
    "\n",
    "values = [(CV1Lo, CV1Hi), (CV2Lo, CV2Hi)]\n",
    "G = [(G11, G12), (G21, G22)]\n",
    "\n",
    "sliders = [widgets.FloatRangeSlider(\n",
    "                value=[CV1Lo, CV1Hi], min=-1.5, max=1.5, step=stepsize,\n",
    "                description='CV1 Limits'),\n",
    "            widgets.FloatRangeSlider(\n",
    "                value=[CV2Lo, CV2Hi], min=-3, max=3.5, step=stepsize,\n",
    "                description='CV2 Limits'),\n",
    "            widgets.FloatSlider(\n",
    "                value=cost_MV1, min=-2, max=2, step=stepsize,\n",
    "                description='MV1 Cost'),\n",
    "            widgets.FloatSlider(\n",
    "                value=cost_MV2, min=-2, max=2, step=stepsize,\n",
    "                description='MV2 Cost'),\n",
    "            widgets.FloatRangeSlider(\n",
    "                value=[MV1_lims[0], MV1_lims[1]], min=-limits, max=limits, step=stepsize,\n",
    "                description='MV1 Limits'),\n",
    "            widgets.FloatRangeSlider(\n",
    "                value=[MV2_lims[0], MV2_lims[1]], min=-limits, max=limits, step=stepsize,\n",
    "                description='MV2 Limits')]\n",
    "\n",
    "output = widgets.Output()\n",
    "with output:\n",
    "    fig = plt.figure(figsize=(6,6))\n",
    "    plt.show(fig)\n",
    "ax = fig.add_subplot(111)\n",
    "ax.grid(True, which='both')\n",
    "plt.xlim((-limits, limits))\n",
    "plt.ylim((-limits, limits))\n",
    "plt.plot(0,0,'kx');\n",
    "# plt.xlabel(r'$\\Delta MV_1$')\n",
    "# plt.ylabel(r'$\\Delta MV_2$')\n",
    "im = plt.imshow((c1 & c2 & c3 & c4).astype(int), extent=(x.min(),x.max(),y.min(),y.max()), origin=\"lower\", cmap=\"Blues\", alpha=0.1);\n",
    "\n",
    "# plot vector field to show the direction of optimization (direction of decreasing cost)\n",
    "q = plt.quiver(xv, yv,-cost_MV1,-cost_MV2, z_obj, cmap='gray', headwidth=4, width=0.003, scale=40, alpha=0.5)\n",
    "\n",
    "nCVs = 2 # How many CV lines?\n",
    "CV_lines_lo = []\n",
    "CV_lines_hi = []\n",
    "\n",
    "colors = ['r', 'b', 'y']\n",
    "for i in range(nCVs):\n",
    "    y_lo = (values[i][0] - G[i][0]*d)/G[i][1]\n",
    "    y_hi = (values[i][1] - G[i][0]*d)/G[i][1]    \n",
    "    line_lo, = ax.plot(d, y_lo, f'--{colors[i]}', label=f'$CV_{i}$ - Lo');\n",
    "    line_hi, = ax.plot(d, y_hi, f'-{colors[i]}', label=f'$CV_{i}$ - Hi');\n",
    "    CV_lines_lo.append(line_lo)\n",
    "    CV_lines_hi.append(line_hi)\n",
    "    \n",
    "mv1_line_lo, = ax.plot(MV1_lims[0] * np.ones(d.shape), d,'g--', label='$MV_1$ - Lo');\n",
    "mv1_line_hi, = ax.plot(MV1_lims[1] * np.ones(d.shape), d,'g', label='$MV_1$ - Hi');\n",
    "mv2_line_lo, = ax.plot(d,MV2_lims[0] * np.ones(d.shape), 'g--', label='$MV_2$ - Lo');\n",
    "mv2_line_hi, = ax.plot(d,MV2_lims[1] * np.ones(d.shape), 'g', label='$MV_2$ - Hi');\n",
    "\n",
    "soln_marker, = ax.plot([0], [0], 'ok')\n",
    "soln_text = ax.text(0, 0, f'(0,0)')\n",
    "\n",
    "obj_func, = ax.plot([], [], '--k')\n",
    "plt.legend(loc='lower center', ncol=4, bbox_to_anchor=(0.5, 1.01))\n",
    "ax.relim()\n",
    "ax.autoscale_view()\n",
    "\n",
    "# Solve the LP\n",
    "def run_lp(values_dict):\n",
    "    '''\n",
    "    Run pulp solver for LP problem with the same gain matrix and CV limits as the example. Returns a tuple with:\n",
    "    (LP solution array [MV1_solution, MV2_solution], objective_value)\n",
    "    '''    \n",
    "    prob = LpProblem(\"DMC_problem\",LpMinimize)\n",
    "    \n",
    "    nMVs = 2\n",
    "    nCVs = 2\n",
    "    \n",
    "    # How many MVs\n",
    "    MVs = []\n",
    "    for i in range(nMVs):\n",
    "        MVs.append(LpVariable(f\"MV{i+1}\",-limits))\n",
    "    \n",
    "    # the objective function\n",
    "    obj = 0\n",
    "    for indx, MV in enumerate(MVs):\n",
    "#         print([f\"MV{indx+1} Cost\"])\n",
    "        obj += values_dict[f\"MV{indx+1} Cost\"]*MV\n",
    "        \n",
    "    prob += obj, \"Cost function of MVs\"\n",
    "    \n",
    "    # constraint formulation in terms of MV1 and MV2\n",
    "    CV_contraint_lo = []\n",
    "    CV_contraint_hi = []\n",
    "    \n",
    "    for i in range(nCVs):\n",
    "        c = 0\n",
    "        for indx, MV in enumerate(MVs):\n",
    "            c += G[i][indx]*MV\n",
    "        prob += c <= values_dict[f'CV{i+1} Limits'][1], f'CV{i+1} High Limit'\n",
    "        prob += c >= values_dict[f'CV{i+1} Limits'][0], f'CV{i+1} Low Limit'\n",
    "        \n",
    "    for indx, MV in enumerate(MVs):\n",
    "        prob += MV <= values_dict[f'MV{indx+1} Limits'][1], f'MV{indx+1} High Limit'\n",
    "        prob += MV >= values_dict[f'MV{indx+1} Limits'][0], f'MV{indx+1} Low Limit'        \n",
    "    \n",
    "    if (prob.solve(PULP_CBC_CMD(msg=0)) == 1):\n",
    "        return [v.varValue for v in prob.variables()], value(prob.objective)\n",
    "    else:\n",
    "        print(\"NOT SOLVED - Infeasibility!\")\n",
    "        return np.zeros(nMVs), 0\n",
    "\n",
    "def handle_slider_change(change):\n",
    "    ## grab slider values for CV contraints\n",
    "    values = [slider.value for slider in sliders]\n",
    "    values_dict = {}\n",
    "    values_dict['CV1 Limits'] = values[0]\n",
    "    values_dict['CV2 Limits'] = values[1]\n",
    "    values_dict['MV1 Limits'] = values[4]\n",
    "    values_dict['MV2 Limits'] = values[5]\n",
    "    values_dict['MV1 Cost'] = values[2]\n",
    "    values_dict['MV2 Cost'] = values[3]\n",
    "    \n",
    "    cost_MV1 = values[2]\n",
    "    cost_MV2 = values[3]\n",
    "    \n",
    "    if (cost_MV2 == 0):\n",
    "        cost_MV2 = 1e-10\n",
    "    \n",
    "    # Update CV constraint lines\n",
    "    for i in range(nCVs):\n",
    "        y_lo = (values[i][0] - G[i][0]*d)/G[i][1]\n",
    "        y_hi = (values[i][1] - G[i][0]*d)/G[i][1]\n",
    "        CV_lines_lo[i].set_data(d, y_lo)\n",
    "        CV_lines_hi[i].set_data(d, y_hi)\n",
    "    \n",
    "    # Recalculate shaded regions\n",
    "    c1 = G11*x+G12*y <= values[0][1]\n",
    "    c2 = G11*x+G12*y >= values[0][0]\n",
    "    c3 = G21*x+G22*y <= values[1][1]\n",
    "    c4 = G21*x+G22*y >= values[1][0]\n",
    "    \n",
    "    # MVs\n",
    "    m1 = (x >= values[4][0]) & (x <= values[4][1])\n",
    "    m2 = (y >= values[5][0]) & (y <= values[5][1])\n",
    "    mv1_line_lo.set_data(values[4][0] * np.ones(d.shape),d);\n",
    "    mv1_line_hi.set_data(values[4][1] * np.ones(d.shape),d);\n",
    "    mv2_line_lo.set_data(d,values[5][0] * np.ones(d.shape));\n",
    "    mv2_line_hi.set_data(d,values[5][1] * np.ones(d.shape));\n",
    "\n",
    "    # Shade the right regions\n",
    "    im.set_data((c1 & c2 & c3 & c4 & m1 & m2).astype(float))\n",
    "    \n",
    "    z_obj = (cost_MV1 * xv) + (cost_MV2 * yv)\n",
    "    q.set_UVC(-cost_MV1, -cost_MV2, z_obj)\n",
    "    \n",
    "    # Find the LP soln\n",
    "    soln, V = run_lp(values_dict)\n",
    "    \n",
    "    if (V == None):\n",
    "        V = 0\n",
    "\n",
    "    y_obj = (1/cost_MV2)*V - (cost_MV1/cost_MV2)*d\n",
    "    v_new = V\n",
    "    \n",
    "    def get_sign(x):\n",
    "        return ['-' if x < 0 else ('-' if x > 0 else '')]\n",
    "\n",
    "    # the soln\n",
    "    soln_marker.set_data(soln[0], soln[1]);\n",
    "    soln_text.set_position((soln[0], soln[1]));\n",
    "    soln_text.set_text(\"({:.1f}, {:.1f})\".format(soln[0], soln[1]))\n",
    "    obj_func.set_data(d, y_obj);\n",
    "    fig.canvas.draw()\n",
    "\n",
    "    with output_text:\n",
    "        clear_output()\n",
    "        print('Cost of current operating point: ${:.2f}'.format(v_new))\n",
    "        \n",
    "    \n",
    "# register slides\n",
    "for widget in sliders:\n",
    "    widget.observe(handle_slider_change, names='value')\n",
    "\n",
    "output_text = widgets.Output()\n",
    "\n",
    "display(widgets.HBox([output, \n",
    "             widgets.VBox([widgets.HTMLMath(\n",
    "                            value=r\"<p>Use the controls to interact with this linear program</p><br><br>\"), \n",
    "                           *sliders,\n",
    "                          output_text])]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "03b0ce1f",
   "metadata": {},
   "source": [
    "## Animated plot\n",
    "<!-- Initialize plot (TODO clean up) -->"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b54187b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Initialize the plot\n",
    "def init_func():\n",
    "    limits = 10\n",
    "    d = np.linspace(-limits, limits, 100)\n",
    "    x,y = np.meshgrid(d,d)\n",
    "    stepsize = 0.1\n",
    "    v_new = 0\n",
    "    \n",
    "    # mesh for vector field\n",
    "    dvec = np.linspace(-limits + (0.1*limits), limits - (0.1*limits), 12)\n",
    "    xv, yv = np.meshgrid(dvec, dvec)\n",
    "    \n",
    "    # constraint formulation in terms of x and y\n",
    "    c1 = G11*x+G12*y <= CV1Hi\n",
    "    c2 = G11*x+G12*y >= CV1Lo\n",
    "    c3 = G21*x+G22*y <= CV2Hi\n",
    "    c4 = G21*x+G22*y >= CV2Lo\n",
    "    \n",
    "    values = [(CV1Lo, CV1Hi), (CV2Lo, CV2Hi)]\n",
    "    G = [(G11, G12), (G21, G22)]\n",
    "    \n",
    "    fig_anim = plt.figure(figsize=(6,6))\n",
    "    ax = fig_anim.add_subplot(111)\n",
    "    ax.grid(True, which='both')\n",
    "    plt.xlim((-limits, limits))\n",
    "    plt.ylim((-limits, limits))\n",
    "    plt.plot(0,0,'kx');\n",
    "    plt.xlabel(r'$\\Delta MV_1$')\n",
    "    plt.ylabel(r'$\\Delta MV_2$')\n",
    "    \n",
    "    im = plt.imshow((c1 & c2 & c3 & c4).astype(int), extent=(x.min(),x.max(),y.min(),y.max()), origin=\"lower\", cmap=\"Blues\", alpha=0.1);\n",
    "    q = plt.quiver(xv, yv,-cost_MV1,-cost_MV2, z_obj, cmap='gray', headwidth=4, width=0.003, scale=40, alpha=0.5)\n",
    " \n",
    "    nCVs = 2 # How many CV lines?\n",
    "    CV_lines_lo = []\n",
    "    CV_lines_hi = []\n",
    "\n",
    "    colors = ['r', 'b', 'y']\n",
    "    for i in range(nCVs):\n",
    "        y_lo = (values[i][0] - G[i][0]*d)/G[i][1]\n",
    "        y_hi = (values[i][1] - G[i][0]*d)/G[i][1]    \n",
    "        line_lo, = ax.plot(d, y_lo, f'--{colors[i]}', label=f'$CV_{i}$ - Lo');\n",
    "        line_hi, = ax.plot(d, y_hi, f'-{colors[i]}', label=f'$CV_{i}$ - Hi');\n",
    "        CV_lines_lo.append(line_lo)\n",
    "        CV_lines_hi.append(line_hi)\n",
    "\n",
    "    soln_marker, = ax.plot([0], [0], 'ok')\n",
    "    soln_text = ax.text(0, 0, f'(0,0)')\n",
    "\n",
    "    obj_func, = ax.plot([], [], '--k')\n",
    "#     plt.legend(loc='lower center', ncol=4, bbox_to_anchor=(0.5, 1.05))\n",
    "    ax.relim()\n",
    "    ax.autoscale_view()\n",
    "    return fig_anim, ax, CV_lines_lo, CV_lines_hi, soln_marker, soln_text, obj_func, im, q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e109db2",
   "metadata": {},
   "source": [
    "## Set up the animation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b79e78d9",
   "metadata": {},
   "outputs": [],
   "source": [
    "n_frames = 50\n",
    "cost_MV1_arr = np.hstack((np.linspace(-1,1,n_frames), np.linspace(1,2,n_frames)))\n",
    "cost_MV2_arr = np.hstack((np.linspace(-1,-1,n_frames), np.linspace(-1,2,n_frames)))\n",
    "n_frames = len(cost_MV1_arr)\n",
    "fig, ax, CV_lines_lo, CV_lines_hi, soln_marker, soln_text, obj_func, im, q = init_func()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "94dd4426",
   "metadata": {},
   "outputs": [],
   "source": [
    "def drawframe(n):\n",
    "    cost_MV1 = cost_MV1_arr[n]\n",
    "    cost_MV2 = cost_MV2_arr[n]\n",
    "    \n",
    "    values = [(CV1Lo, CV1Hi), (CV2Lo, CV2Hi), cost_MV1, cost_MV2]\n",
    "    values_dict = {}\n",
    "    values_dict['CV1 Limits'] = (CV1Lo, CV1Hi)\n",
    "    values_dict['CV2 Limits'] = (CV2Lo, CV2Hi)\n",
    "    values_dict['MV1 Limits'] = (-limits, limits)\n",
    "    values_dict['MV2 Limits'] = (-limits, limits)\n",
    "    values_dict['MV1 Cost'] = cost_MV1\n",
    "    values_dict['MV2 Cost'] = cost_MV2\n",
    "    \n",
    "    if (cost_MV2 == 0):\n",
    "        cost_MV2 = 1e-10\n",
    "    \n",
    "    # Update CV constraint lines\n",
    "    for i in range(nCVs):\n",
    "        y_lo = (values[i][0] - G[i][0]*d)/G[i][1]\n",
    "        y_hi = (values[i][1] - G[i][0]*d)/G[i][1]\n",
    "        CV_lines_lo[i].set_data(d, y_lo)\n",
    "        CV_lines_hi[i].set_data(d, y_hi)\n",
    "    \n",
    "    # Recalculate shaded regions\n",
    "    c1 = G11*x+G12*y <= values[0][1]\n",
    "    c2 = G11*x+G12*y >= values[0][0]\n",
    "    c3 = G21*x+G22*y <= values[1][1]\n",
    "    c4 = G21*x+G22*y >= values[1][0]\n",
    "\n",
    "    # Shade the right regions\n",
    "    im.set_data((c1 & c2 & c3 & c4).astype(float))\n",
    "    \n",
    "    z_obj = (cost_MV1 * xv) + (cost_MV2 * yv)\n",
    "    q.set_UVC(-cost_MV1, -cost_MV2, z_obj)\n",
    "    \n",
    "    # Find the LP soln\n",
    "    soln, V = run_lp(values_dict)    \n",
    "    if (V == None):\n",
    "        V = 0\n",
    "\n",
    "    y_obj = (1/cost_MV2)*V - (cost_MV1/cost_MV2)*d\n",
    "    v_new = V\n",
    "    \n",
    "    # the soln\n",
    "    soln_marker.set_data(soln[0], soln[1]);\n",
    "    soln_text.set_position((soln[0], soln[1]));\n",
    "    soln_text.set_text(\"({:.1f}, {:.1f})\".format(soln[0], soln[1]))\n",
    "    obj_func.set_data(d, y_obj);\n",
    "    \n",
    "    ax.set_title('Cost MV1 = {:.2f}, Cost MV2 = {:.2f}, Cost Function = ${:.2f}'.format(cost_MV1,cost_MV2, V), fontsize=10)\n",
    "    \n",
    "    # Print out the current frame using \"\\r\", the 'carriage return' character, as our end character.\n",
    "    # This makes Python print the frame on the same line.\n",
    "    print(\"Frame: {:0d}\".format(n), end=\"\\r\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dea575d4",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import animation\n",
    "anim = animation.FuncAnimation(fig, drawframe, frames=n_frames, interval=100)\n",
    "from IPython.display import HTML\n",
    "HTML(anim.to_jshtml())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "defa430a",
   "metadata": {},
   "source": [
    "## A note on the role of MV vs. CV constraints in the LP solution: \n",
    "\n",
    "As seen in the graphs above, the LP solution will be at the intersection of active constraints that results in the lowest cost, and this can be a combination of both CV and MV constraints. \n",
    "\n",
    "> If the LP solution lies at both CV constraints (i.e. both CVs are constrained), then the MVs have freedom to move as their limits are \"wider\"\n",
    "\n",
    "> If the LP solution lies at both MV constraints, then the MVs cannot move and the CV constraints are not active\n",
    "\n",
    "We can generalize this idea and say that:\n",
    "\n",
    "#####  <center>\\# (Active CV constraints)  = \\# (Available MVs) = \\# (Degrees of Freedom in the system)</center>\n",
    "\n",
    "So:\n",
    "- if both CVs are constrained, then we have 2 active CV constraints and 2 available MVs &rarr; 2 Degrees of Freedom\n",
    "- if both MVs are constrained then we have 0 active CV constraints and 0 available MVs &rarr; 0 Degrees of Freedom\n",
    "\n",
    "These definitions tie into how we configure gain matrices in practice. We use the set of \"active CV constraints\" and \"available MVs\" as the variables that we \"pivot\" the gain matrix around when analyzing the LP. \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}