157_2016_Проектна задача по Операциони истражувања.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Проектна задача по предметот Операциони истражувања\n",
    "\n",
    "## Благој Христов 157/2016\n",
    "\n",
    "\n",
    "## Нелинеарно програмирање во *Python* - решение на два примери користејќи ја библиотеката *SciPy*"
   ]
  },
  {
   "attachments": {
    "pr1.jpg": {
     "image/jpeg": "/9j/4AAQSkZJRgABAQEAYABgAAD/2wBDAAMCAgMCAgMDAwMEAwMEBQgFBQQEBQoHBwYIDAoMDAsKCwsNDhIQDQ4RDgsLEBYQERMUFRUVDA8XGBYUGBIUFRT/2wBDAQMEBAUEBQkFBQkUDQsNFBQUFBQUFBQUFBQUFBQUFBQUFBQUFBQUFBQUFBQUFBQUFBQUFBQUFBQUFBQUFBQUFBT/wAARCAEJAYwDASIAAhEBAxEB/8QAHwAAAQUBAQEBAQEAAAAAAAAAAAECAwQFBgcICQoL/8QAtRAAAgEDAwIEAwUFBAQAAAF9AQIDAAQRBRIhMUEGE1FhByJxFDKBkaEII0KxwRVS0fAkM2JyggkKFhcYGRolJicoKSo0NTY3ODk6Q0RFRkdISUpTVFVWV1hZWmNkZWZnaGlqc3R1dnd4eXqDhIWGh4iJipKTlJWWl5iZmqKjpKWmp6ipqrKztLW2t7i5usLDxMXGx8jJytLT1NXW19jZ2uHi4+Tl5ufo6erx8vP09fb3+Pn6/8QAHwEAAwEBAQEBAQEBAQAAAAAAAAECAwQFBgcICQoL/8QAtREAAgECBAQDBAcFBAQAAQJ3AAECAxEEBSExBhJBUQdhcRMiMoEIFEKRobHBCSMzUvAVYnLRChYkNOEl8RcYGRomJygpKjU2Nzg5OkNERUZHSElKU1RVVldYWVpjZGVmZ2hpanN0dXZ3eHl6goOEhYaHiImKkpOUlZaXmJmaoqOkpaanqKmqsrO0tba3uLm6wsPExcbHyMnK0tPU1dbX2Nna4uPk5ebn6Onq8vP09fb3+Pn6/9oADAMBAAIRAxEAPwD7p+Pnxk8YfD/xp8O/Cfgjw3pfiHW/F0t/HGNXvntIYRbQCZiWVGPK7u3UCsb/AITD9pn/AKJv4B/8Ke4/+RqT43f8nXfs2f8AXx4i/wDTYa+hKAPnz/hMP2mf+ib+Af8Awp7j/wCRqP8AhMP2mf8Aom/gH/wp7j/5Gr6DooA+fP8AhMP2mf8Aom/gH/wp7j/5Go/4TD9pn/om/gH/AMKe4/8AkavoOigD58/4TD9pn/om/gH/AMKe4/8Akaj/AITD9pn/AKJv4B/8Ke4/+Rq+g6KAPnz/AITD9pn/AKJv4B/8Ke4/+RqP+Ew/aZ/6Jv4B/wDCnuP/AJGr6DooA+fP+Ew/aZ/6Jv4B/wDCnuP/AJGo/wCEw/aZ/wCib+Af/CnuP/kavoOigD58/wCEw/aZ/wCib+Af/CnuP/kaj/hMP2mf+ib+Af8Awp7j/wCRq+g6KAPnz/hMP2mf+ib+Af8Awp7j/wCRqP8AhMP2mf8Aom/gH/wp7j/5Gr6DooA+fP8AhMP2mf8Aom/gH/wp7j/5Go/4TD9pn/om/gH/AMKe4/8AkavoOigD40+Mv7Wnxj+Aq+Gj4x8FeAdPHiDUU02z2+JLhsMcbpH/ANH+WNMrubtuHrXpH/CYftM/9E38A/8AhT3H/wAjVxv7b37MHhH42a/8OtY8R3WsGY67ZaH5FrdiOFbaZ3aXClDh2IXLDnCqO1fUfhvQ4vDHh3TNHguLm7g0+2jtY57yTzJnVFCgu2BubAGT3PNAHh3/AAmH7TP/AETfwD/4U9x/8jUf8Jh+0z/0TfwD/wCFPcf/ACNX0HRQB8+f8Jh+0z/0TfwD/wCFPcf/ACNR/wAJh+0z/wBE38A/+FPcf/I1fQdFAHz5/wAJh+0z/wBE38A/+FPcf/I1H/CYftM/9E38A/8AhT3H/wAjV9B0UAfPn/CYftM/9E38A/8AhT3H/wAjUf8ACYftM/8ARN/AP/hT3H/yNX0HRQB8+f8ACYftM/8ARN/AP/hT3H/yNR/wmH7TP/RN/AP/AIU9x/8AI1fQdFAHz5/wmH7TP/RN/AP/AIU9x/8AI1H/AAmH7TP/AETfwD/4U9x/8jV9B0UAfPn/AAmH7TP/AETfwD/4U9x/8jUf8Jh+0z/0TfwD/wCFPcf/ACNX0HRQB8+f8Jh+0z/0TfwD/wCFPcf/ACNR/wAJh+0z/wBE38A/+FPcf/I1fQdFAHz5/wAJh+0z/wBE38A/+FPcf/I1H/CYftM/9E38A/8AhT3H/wAjV9B0UAfPn/CYftM/9E38A/8AhT3H/wAjUf8ACYftM/8ARN/AP/hT3H/yNX0HRQB8+f8ACYftM/8ARN/AP/hT3H/yNR/wmH7TP/RN/AP/AIU9x/8AI1fQdFAHz5/wmH7TP/RN/AP/AIU9x/8AI1H/AAmH7TP/AETfwD/4U9x/8jV9B0UAfPn/AAmH7TP/AETfwD/4U9x/8jUf8Jh+0z/0TfwD/wCFPcf/ACNX0HRQB8+f8Jh+0z/0TfwD/wCFPcf/ACNR/wAJh+0z/wBE38A/+FPcf/I1fQdFAHzT4o+KH7SHhPwzq2uXnw08CyWmmWk17MkPiacyMkaF2Cg24BOFOMkfWva/hN40l+JHwr8G+LZrVLKbXtFstVe1jYssLTwJKUBPUAvjPtVf41f8kb8ef9gC/wD/AEnkrD/Zc/5Nl+Ef/YoaR/6RRUAcV8bv+Trv2bP+vjxF/wCmw19CV89/G7/k679mz/r48Rf+mw19CUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQB5P+0F934cf9jrpn85K9Yryf9oL7vw4/wCx10z+clesUAFFeU/GT4tX3wn8WeB5rlLUeENVmvLPU7iVG82CZLZ54CrbsBT5UikEHquCK8o+D37UnjP4hR+CtI1PS9K03xVf6zdjWLZY5BHb6XFaLdpIoL5V2S4tV3MSMs3FAH1bUa3EUm3bIjbslcMDnHXFfJPgj9qbxPqXxJ8D2N1rvh/xHp3ijU5tNubHQ9DvFt9MIilkiMeps5huW/dYYBRnJKgAV47oPxgsvgX4U+Bni29sLrVPs2i+LEhs7SNmMkhvkI3MAfLjABZnPCqpPPQgH6OUV8w+Lv2hPHHwrih/4SZvDurPqfg661nTZtIt5Vgk1OOVAlsrNIxkjKzw4PBbDHjIAZ4J/ad8V+KLDTVns9BstT0vwvqWq+KFupGgtra+t7hraKISsx8qNpYp2Jbcdq0AfUNFfHvhX9pX4g654im8O6VrOg+MdT1Pwxfaxp11B4ZvtOtLe+txGRbpJLJi7ibzMb0II4J+8BXpvw3/AGgb/wCLPjTwJZaDBaDSL7wsfEWvPIjNJbPIwigt0O7Ct5qz53A8RdutAHutFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFAHGfGr/kjfjz/ALAF/wD+k8lYf7Ln/Jsvwj/7FDSP/SKKtz41f8kb8ef9gC//APSeSsP9lz/k2X4R/wDYoaR/6RRUAcV8bv8Ak679mz/r48Rf+mw19CV89/G7/k679mz/AK+PEX/psNfQlABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAeT/ALQX3fhx/wBjrpn85K9Yr5P/AG2P2mvBnwZ174daP4lOpQ3S6/Y62Xt7NpIvssTyLKwfOCynGUHzfMpxyK+nvDevW/ijw9pms2kc8VrqFtHdwpdRGKUI6hl3oeVbBGQeR0NAHLfGr4N6D8ePAdx4T8RSXkGnzTRXHn6fKsc8bxsGBVmVgMjKnjoxFZll+zv4RsfibrnjiOO7Op6voyaHPatKv2VIAqqWRQoZXZUjUtu5Ea8DFenUUAeA6J+yHY6Ovg2NviD4vvbXwbeR3Hh+1nks/Kso1Vk8ohbceaDGxTdJucD7rKck9L4V/Zn8J+FV8KJFPqN/B4dstS0+3gvpInjuIb6USTiYCMbuRhcbRgnIbrXrNFAHiL/sj+Drjwp4I8O3Wo65ead4N1catpHn3UbPGoYstozeX81uPlAU/NhFG7irn/DK3glofijE/wDaMkfxEOdVU3CjyBhuLfCjYN7vJg7vmYnpxXsVFAHiGm/s9t4S8UaN47vfGvjHxl4h8OWNxa28M5sR9qtXQZthEsMaBiY0O/KuWA3PjgVf2T/g9N8ONJ8X6/qGjz6BqXivWrjUY9Iupkll06z8xzBbkozKCN8jlVJAMpGTiveaKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigDjPjV/wAkb8ef9gC//wDSeSsP9lz/AJNl+Ef/AGKGkf8ApFFW58av+SN+PP8AsAX/AP6TyVh/suf8my/CP/sUNI/9IoqAOK+N3/J137Nn/Xx4i/8ATYa+hK+e/jd/ydd+zZ/18eIv/TYa+hKACiiigAooooAKKKKACiiigAooooA8j/aQ+Mmu/BvQ/CMnhrw3beKda8SeI7Xw9bWV3ffY41eaOZxI0mxuAYsYx0J9MHmP+Fh/tIf9Ee8Jf+Fi3/yNR+13/wAhD4D/APZT9K/9J7uvoOgD58/4WH+0h/0R7wl/4WLf/I1H/Cw/2kP+iPeEv/Cxb/5Gr6DooA+fP+Fh/tIf9Ee8Jf8AhYt/8jUf8LD/AGkP+iPeEv8AwsW/+Rq+g6KAPiT4+fC34x/tFDwiPFPwX8Jyf8I5q0epw7PF+fOQf6y2fNr/AKuTCbgOfkFer/8ACwv2kBwPg94R/wDCxb/5Gr6EooA+fP8AhYf7SH/RHvCX/hYt/wDI1H/Cw/2kP+iPeEv/AAsW/wDkavoOigD58/4WH+0h/wBEe8Jf+Fi3/wAjUf8ACw/2kP8Aoj3hL/wsW/8AkavoOigD58/4WH+0h/0R7wl/4WLf/I1H/Cw/2kP+iPeEv/Cxb/5Gr6DooA+fP+Fh/tIf9Ee8Jf8AhYt/8jUf8LD/AGkP+iPeEv8AwsW/+Rq+g6KAPnz/AIWH+0h/0R7wl/4WLf8AyNR/wsP9pD/oj3hL/wALFv8A5Gr6DooA+fP+Fh/tIf8ARHvCX/hYt/8AI1H/AAsP9pD/AKI94S/8LFv/AJGr6DooA+fP+Fh/tIf9Ee8Jf+Fi3/yNR/wsP9pD/oj3hL/wsW/+Rq+g6KAPnz/hYf7SH/RHvCX/AIWLf/I1H/Cw/wBpD/oj3hL/AMLFv/kavoOigD5W+Jf7Qvx7+FHgPWvF2ufB3wydI0iD7TdfZfFpeQRggEhfs4zjOcZ7V9P6XeHUNNtLopsM8KS7c5xuUHH614v+3B/yaT8Uv+wLL/Na9g8M/wDIt6V/16Rf+gCgDJ8a6p4x037H/wAIp4d0jXg+/wC0f2prUmn+VjG3bstpt+ctnO3GB1zx51oHxY+KniLXPE2lW3w48MJceH72Oxumk8YzhWd7aG4BTGncjZOg5xyD9a9trzD4V/8AJTfjP/2MVn/6ZtPoA7Hwbf8AiXUNNlk8UaLpuh34lKx2+mam9/G0eBhi7wQkHO4bdp4AOecDeoooAKKKKACiiigAooooA4z41f8AJG/Hn/YAv/8A0nkrD/Zc/wCTZfhH/wBihpH/AKRRVufGr/kjfjz/ALAF/wD+k8lYf7Ln/Jsvwj/7FDSP/SKKgDivjd/ydd+zZ/18eIv/AE2GvoSvnv43f8nXfs2f9fHiL/02GvoSgAooooAKKKKACiiigAooooAKKKKAPnz9rv8A5CHwH/7KfpX/AKT3dfQdfPn7Xf8AyEPgP/2U/Sv/AEnu6+g6ACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooA8N/bg/5NJ+KX/YFl/mteweGf+Rb0r/r0i/8AQBXj/wC3B/yaT8Uv+wLL/Na9g8M/8i3pX/XpF/6AKANKvMPhX/yU34z/APYxWf8A6ZtPr0+vMPhX/wAlN+M//YxWf/pm0+gD0+iiigAooooAKKKKACiiigDjPjV/yRvx5/2AL/8A9J5Kw/2XP+TZfhH/ANihpH/pFFW58av+SN+PP+wBf/8ApPJWH+y5/wAmy/CP/sUNI/8ASKKgDivjd/ydd+zZ/wBfHiL/ANNhr6Er57+N3/J137Nn/Xx4i/8ATYa+hKACiiigAooooAKKKKACiiigAooooA+fP2u/+Qh8B/8Asp+lf+k93X0HXz5+13/yEPgP/wBlP0r/ANJ7uvoOgAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKAPDf24P8Ak0n4pf8AYFl/mteweGf+Rb0r/r0i/wDQBXj/AO3B/wAmk/FL/sCy/wA1r2Dwz/yLelf9ekX/AKAKANKvMPhX/wAlN+M//YxWf/pm0+vT68w+Ff8AyU34z/8AYxWf/pm0+gD0+iiigAooooAKKKKACiiigDjPjV/yRvx5/wBgC/8A/SeSsP8AZc/5Nl+Ef/YoaR/6RRVufGr/AJI348/7AF//AOk8lYf7Ln/Jsvwj/wCxQ0j/ANIoqAOK+N3/ACdd+zZ/18eIv/TYa+hK+e/jd/ydd+zZ/wBfHiL/ANNhr6EoAKKKKACiiigAooooAKKKKACiiigD58/a7/5CHwH/AOyn6V/6T3dfQdfPn7Xf/IQ+A/8A2U/Sv/Se7r6DoAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigDw39uD/k0n4pf9gWX+a17B4Z/wCRb0r/AK9Iv/QBXj/7cH/JpPxS/wCwLL/Na9g8M/8AIt6V/wBekX/oAoA0q8w+Ff8AyU34z/8AYxWf/pm0+vT68w+Ff/JTfjP/ANjFZ/8Apm0+gDqfiND4pn8J3q+DtQ0vS9cxmO61e0kuYUXB3Hy0kQlvTJxnqDXzH4f+J3jC4/Y4+E7DxFeDxX4y1Wz0GfxC7b7iET3MglmUtnDiNGCk5wSPSvr64i8+3lizt3qVzjpkYrxG3/ZjS3/Z58PfDZPEkian4flhvdM8Qx2YUw3kM5mil8kuQQCdpXdypIyM8AFf4VTav8Pfj54i+Glx4j1fxPoUmgW/iHT59eujdXdq5naCaIzEbmQkIwDZwd2OKzviX+0t4v8AC/i34kaZ4d8CWGsad4CsLbVNT1C+1lrYzwSQGZo4UWB8ygK2NxC/LyecV2/wt+E2u+HvGmu+N/GniCz8ReLtVtYNOV9MsGs7O0tIizLFHG0kjEs7s7Mz8nAAAFV9f+AX9uXnxgn/ALd8j/hYOmwadt+x7vsHl2rwb/8AWDzc792PkxjGe9AGZ8N/j14k8TfELRPD3ibwfa+HrXxLocmvaJNa6mbuXyUaMNFcp5SLHJtmRsIzqMkZOM17fXmWn/Bf7B468B+I/wC2PM/4Rbw/PoX2b7Lj7V5gg/e7t/yY8j7uGzu68c+m0AFFFFAHGfGr/kjfjz/sAX//AKTyVh/suf8AJsvwj/7FDSP/AEiirc+NX/JG/Hn/AGAL/wD9J5Kw/wBlz/k2X4R/9ihpH/pFFQBxXxu/5Ou/Zs/6+PEX/psNfQlfPfxu/wCTrv2bP+vjxF/6bDX0JQAUUUUAFFFFABRRRQAUUUUAFFFFAHz5+13/AMhD4D/9lP0r/wBJ7uvoOvnz9rv/AJCHwH/7KfpX/pPd19B0AFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQB4b+3B/yaT8Uv+wLL/Na9g8M/8i3pX/XpF/6AK8f/AG4P+TSfil/2BZf5rXsHhn/kW9K/69Iv/QBQBpV5h8K/+Sm/Gf8A7GKz/wDTNp9en15h8K/+Sm/Gf/sYrP8A9M2n0Aen0UUUAFFFFABRRRQAUUUUAcZ8av8Akjfjz/sAX/8A6TyVh/suf8my/CP/ALFDSP8A0iirc+NX/JG/Hn/YAv8A/wBJ5Kw/2XP+TZfhH/2KGkf+kUVAHFfG7/k679mz/r48Rf8ApsNfQlfPfxu/5Ou/Zs/6+PEX/psNfQlABRRRQAUUUUAFFFFABRRRQAUUUUAfPn7Xf/IQ+A//AGU/Sv8A0nu6+g6+fP2u/wDkIfAf/sp+lf8ApPd19B0AFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQB4b+3B/wAmk/FL/sCy/wA1r2Dwz/yLelf9ekX/AKAK8f8A24P+TSfil/2BZf5rXsHhn/kW9K/69Iv/AEAUAaVeYfCv/kpvxn/7GKz/APTNp9en15h8K/8Akpvxn/7GKz/9M2n0Aen0UUUAFFFFABRRRQAUUUUAcZ8av+SN+PP+wBf/APpPJWH+y5/ybL8I/wDsUNI/9Ioq3PjV/wAkb8ef9gC//wDSeSsP9lz/AJNl+Ef/AGKGkf8ApFFQBxXxu/5Ou/Zs/wCvjxF/6bDX0JXz38bv+Trv2bP+vjxF/wCmw19CUAFFFFABRRRQAUUUUAFFFFABRRRQB4H+194c8Wazo3w01Twf4XuPF994b8b2GuXOmWtzFBI1tFDcq5DSsq5zIoxnv9SIP+Ghvin/ANG4+Kv/AAd6Z/8AHq9i8XeM7fwfP4fiuLeWc6zqkelxGPH7t3SRwzZ7fuyOPWuhoA+fP+Ghvin/ANG4+Kv/AAd6Z/8AHqP+Ghvin/0bj4q/8Hemf/Hq+g6KAPnz/hob4p/9G4+Kv/B3pn/x6j/hob4p/wDRuPir/wAHemf/AB6voOigD58/4aG+Kf8A0bj4q/8AB3pn/wAeo/4aG+Kf/RuPir/wd6Z/8er6DooA+fP+Ghvin/0bj4q/8Hemf/HqP+Ghvin/ANG4+Kv/AAd6Z/8AHq+g6KAPnz/hob4p/wDRuPir/wAHemf/AB6j/hob4p/9G4+Kv/B3pn/x6voOigD58/4aG+Kf/RuPir/wd6Z/8eo/4aG+Kf8A0bj4q/8AB3pn/wAer6DooA+fP+Ghvin/ANG4+Kv/AAd6Z/8AHqP+Ghvin/0bj4q/8Hemf/Hq+g6KAPnz/hob4p/9G4+Kv/B3pn/x6j/hob4p/wDRuPir/wAHemf/AB6voOigD58/4aG+Kf8A0bj4q/8AB3pn/wAeo/4aG+Kf/RuPir/wd6Z/8er6DooA+fP+Ghvin/0bj4q/8Hemf/HqP+Ghvin/ANG4+Kv/AAd6Z/8AHq+g6KAPnz/hob4p/wDRuPir/wAHemf/AB6j/hob4p/9G4+Kv/B3pn/x6voOigD46/aC8e/GD4w/Bfxd4L0/9nvxHY3mt2TWcdzc61pxji3EZYhZSTgA8d/brX1voVvJaaHp8Eq7JYreNHX0YKARWZ4+8Z2/w/8ACt1rt3by3MFvJBG0UONx82ZIgRnjguD+FdDQBkeKtavtA0aS807QL7xLdKyhdP02W3jmcE4JDTyxx4A5OWz6ZrxbwX4m8eeHPGHxA1W5+DXip7fxBqtvfWyx6lopZESwtbYh834wd8DnjPBHfgfQFFAHH+DfGuu+JdQmt9U8Aa94Shji8xbrVbrT5Y5GyBsAtrqVg2CTyoGAec4B7CiigAooooAKKKKACiiigDjPjV/yRvx5/wBgC/8A/SeSsP8AZc/5Nl+Ef/YoaR/6RRVufGr/AJI348/7AF//AOk8lYf7Ln/Jsvwj/wCxQ0j/ANIoqAOK+N3/ACdd+zZ/18eIv/TYa+hK+e/jd/ydd+zZ/wBfHiL/ANNhr6EoAKKKKACiiigAooooAKKKKACiiigDx79pBdXaD4dDQnso9U/4S618ltRR2gB+z3P3ghDfka0vI+M//P74E/8AAO9/+O11PjTwUnjK48Nyvdta/wBi6tFqoCx7vNKRyJsPIxnzM556dK6agDzDyPjP/wA/vgT/AMA73/47R5Hxn/5/fAn/AIB3v/x2vT6KAPMPI+M//P74E/8AAO9/+O0eR8Z/+f3wJ/4B3v8A8dr0+igDzDyPjP8A8/vgT/wDvf8A47R5Hxn/AOf3wJ/4B3v/AMdr0+igDzDyPjP/AM/vgT/wDvf/AI7R5Hxn/wCf3wJ/4B3v/wAdr0+igDzDyPjP/wA/vgT/AMA73/47R5Hxn/5/fAn/AIB3v/x2vT6KAPMPI+M//P74E/8AAO9/+O0eR8Z/+f3wJ/4B3v8A8drh/wBobW/iB4f+KnwlbT/EVnp/gvUPFVpYXOn2lvKl7cs0UrMss3mbTF8v3AgycZJxXn/7THxH1HSfjJrmk3us+PF0iw8LpqOm2fw4DyS211vk3y34iBZVIEezzf3ZAfPegD3jyPjP/wA/vgT/AMA73/47R5Hxn/5/fAn/AIB3v/x2pvhT401TUP2ffDninWruz1rWG0CO/u7ixkUw3Eoh3MVZfl5IOcDGc4rx/wAO/tW/EjxJD8PJ4fhfpMMXxCtJJdBSXxEweCSOJZXa7xbHZGULMpj3scAFQW4APWvI+M//AD++BP8AwDvf/jtHkfGf/n98Cf8AgHe//HaufA74oXfxW8JX17qekx6JrWl6pd6NqNlBcG4hW4t5NjGOQqpZGGGBKg84I4r0OgDzDyPjP/z++BP/AADvf/jtHkfGf/n98Cf+Ad7/APHa9PooA8w8j4z/APP74E/8A73/AOO0eR8Z/wDn98Cf+Ad7/wDHa9PooA8w8j4z/wDP74E/8A73/wCO0eR8Z/8An98Cf+Ad7/8AHa9PooA+bPj9D8WB8K9TOqXng1rH7VY+YLS1uxL/AMfkOMFpCOuPwzX0nXM/EbwUnxC8H3egvdtYrcS28hnWPeV8qdJcYyOvl46966agAooooAKKKKACiiigAooooAKKKKAOM+NX/JG/Hn/YAv8A/wBJ5Kw/2XP+TZfhH/2KGkf+kUVbnxq/5I348/7AF/8A+k8lYf7Ln/Jsvwj/AOxQ0j/0iioA4r43f8nXfs2f9fHiL/02GvoSvnv43f8AJ137Nn/Xx4i/9Nhr6EoAKKKKACiiigAooooAK8o/aZm+Jun/AAvn1X4UTWr+JtKuYr99NubYTf2nbRndLapn7rOOhHJxtBUsCPV6KAPOfgD8c/D/AO0N8NdO8XeH3MazZhvbCU/vrG6UDzIJB2ZSeDjkFSODXo1fGXxw8L6t+x/8V7n47eB7Ca98Ca1Ikfj/AMNWa/dBbjUoU6b1LEt7kk8OzL9ceFfFGleNvDem6/oV9Dqej6lAlza3cByksbDII/wPIPBoA1aKKKACiiigAooooAKKKKACiiigAooooA4P4ofC7/hZOoeCLr+0/wCzv+EZ1+DXNv2fzftPlxyJ5Wdy7M+ZndzjHTmuS8WfBfxhb/EjxF4v8A+MNN0GfxNaW9nq1prejtfxgwqyRzwFJ4ijhWI2tuU9cV7TRQBwngX4S2Xw/wDg7Y/D7T7yaW1tNOewW9nUF3Zw26QqMDlnZto45xXMeG/2e/8AhH4Pg1H/AG/9o/4V3aS2ufse3+0N9oLfd/rD5WMbsfP6e9exUUAcP8Kfhn/wrG18TQ/2l/aX9ta/e65u8jyvJ+0MG8r7zbtuPvcZ9BXcUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAcZ8av+SN+PP+wBf/APpPJWH+y5/ybL8I/wDsUNI/9Ioq3PjV/wAkb8ef9gC//wDSeSsP9lz/AJNl+Ef/AGKGkf8ApFFQB51+1Bfar4T+NXwN8Z2nhHxL4t0rQbjWTfxeGdMe+niE9kIYyVXAALN3I4B9K0P+GxIv+iL/ABi/8I+T/wCOV9CUUAfPf/DYkX/RF/jF/wCEfJ/8co/4bEi/6Iv8Yv8Awj5P/jlfQlFAHz3/AMNiRf8ARF/jF/4R8n/xyj/hsSL/AKIv8Yv/AAj5P/jlfQlFAHz3/wANiRf9EX+MX/hHyf8Axyj/AIbEi/6Iv8Yv/CPk/wDjlfQlFAHz3/w2JF/0Rf4xf+EfJ/8AHKP+GxIv+iL/ABi/8I+T/wCOV9CUUAfOt9+1tZanZXFnefBD4vXVpcRtFNBN4Md0kRhhlZS+CCCQQfWvlr4WftFS/sX+Nta0u+8B/ESw+BGt3Rm0SHXdCe2udHvpMs1tEZHCvE2GIG/dxnBO8t+l9cp8Uvhl4f8AjF4D1fwh4oslvtG1OExSp0ZD1WRD/C6sAynsQKAPJh+2LCeR8F/jER/2J8n/AMco/wCGxIv+iL/GL/wj5P8A45XGfsw/E7xB8HvHsn7O3xSvWudYsYjL4P8AEc3Ca5py52xkn/ltGoIxnOFI52hn+t6APnv/AIbEi/6Iv8Yv/CPk/wDjlH/DYkX/AERf4xf+EfJ/8cr6EooA+e/+GxIv+iL/ABi/8I+T/wCOUf8ADYkX/RF/jF/4R8n/AMcr6EooA+e/+GxIv+iL/GL/AMI+T/45R/w2JF/0Rf4xf+EfJ/8AHK+hKKAPnv8A4bEi/wCiL/GL/wAI+T/45R/w2JF/0Rf4xf8AhHyf/HK+hKKAPnv/AIbEi/6Iv8Yv/CPk/wDjlH/DYkX/AERf4xf+EfJ/8cr6EooA+e/+GxIv+iL/ABi/8I+T/wCOUf8ADYkX/RF/jF/4R8n/AMcr6EooA+e/+GxIv+iL/GL/AMI+T/45R/w2JF/0Rf4xf+EfJ/8AHK+hKKAPnv8A4bEi/wCiL/GL/wAI+T/45R/w2JF/0Rf4xf8AhHyf/HK+hKKAPnv/AIbEi/6Iv8Yv/CPk/wDjlH/DYkX/AERf4xf+EfJ/8cr6EooA+e/+GxIv+iL/ABi/8I+T/wCOUf8ADYkX/RF/jF/4R8n/AMcr6EooA+e/+GxIv+iL/GL/AMI+T/45R/w2JF/0Rf4xf+EfJ/8AHK+hKKAPnv8A4bEi/wCiL/GL/wAI+T/45R/w2JF/0Rf4xf8AhHyf/HK+hKKAPnv/AIbEi/6Iv8Yv/CPk/wDjlH/DYkX/AERf4xf+EfJ/8cr6EooA+e/+GxIv+iL/ABi/8I+T/wCOUf8ADYkX/RF/jF/4R8n/AMcr6EooA+e/+GxIv+iL/GL/AMI+T/45R/w2JF/0Rf4xf+EfJ/8AHK+hKKAPnv8A4bEi/wCiL/GL/wAI+T/45R/w2JF/0Rf4xf8AhHyf/HK+hKKAPnv/AIbEi/6Iv8Yv/CPk/wDjlH/DYkX/AERf4xf+EfJ/8cr6EooA+WfiP+1NJ4p+HvijRbL4L/F/7ZqOl3VnB5nhCRV8ySFkXJ38DJFexfs56Te6D+z38MNM1K0m0/UbLwtpdtc2lzGY5YJUtIldHU8qysCCDyCK9EooAKKKKACiiigAooooAKKKKACiiigAooooA8d/ae/Z3sf2hvAcdlFdHRPF2kTDUPDviCElZtPvFwVYMOdjFVDAegI5VSML9k/9oi++LGk6r4Q8bWg0P4s+EXFn4g0lwF83HC3cQHBjkGDleATx8rIT7/XzL+1h8C/EF5q2lfGX4Vqtt8VvCiE/ZlHya7YjmSzlA+8cZ2d+ccHaygH01RXnH7P/AMdPD/7Q/wANdP8AF2gMYvMzBfafKf31hdKB5kEg9VJ4OBlSp716PQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQB8Z/HPwrq37IfxWufjx4GsJr3wPrDrH8QPDVoP4S3GpQp03qWJb3JJ4d2X628KeKtJ8ceG9N8QaFfQ6no+pQLc2t3AcpLGwyCPT3B5ByDzWheWdvqVnPaXcEdza3EbRSwTIHSRGGGVlPBBBIIPrXxX4bvLj9gD4wReFtSmlf4A+NL5jot/OxZPDWoOSWtpGP3YXOSCen3v4ZGIB9tUUgIYAg5FLQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABRRRQAUUUUAFFFFABXLfE/4aeH/AIweBNX8IeKLJb/RdUhMU0fRkPVZEP8AC6sAyt2IFdTUdxbpdW8sEo3RSKUZckZBGDyKAPzCj/bo8W/sqw6x8CLprPx14n0DVIdJ0XxV54eAWL42idActPEGVNucA5DH938/6g18v+I/+CbXwH1zWINVtPC9z4e1CGdbgTaTqEyKXVgwzHIXQDI7KK+oKACoba8gvEZ7eaOdVdo2aNwwDA4KnHcHqKfKhkidFdoywIDrjK+4yMZr4G8I+L/Fvws+HI0jQ/EHiLUr3xX8StX0p7u2srK4vLWKKa5kle2jMaRmaby+TJuVSSVUYC0AfflFfFGofFj4saPoN/oFvrGt6ZcL4s0TT9L1zxZp9kdR+y3jMssdxDD+7baVyr7ULBh6E13K6l430r426v4cvfifew+FPB/hmx12+ubyxtTLfymW5LmZ1iGyIrGN6xgHCgKV5yAfT1FfFvw9+MvjjT/iJ8NLi61nxl4g0TxXHe/bpvEGl2Njp10UtHuIn0+JALmNfk48wEFWGSTXT/DXxV491zwT4L+Luo/E2FLfXpHnufB99DaxWD27LIY7SzYRiY3K7V5Z23kMCAKAPquiviT4U/Gn4o+Im+HfjOeXxdfW/ibUIE1LTb/TdPt9Bis5yQv2SRT5++PKkFmYybWyBxj7boAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAK8yvP2bvh/fadrtjLo9wLbWNUGtTrHqd0hgvslvtFsRKDbSFmJJh2ZJ5zXptFAHm2j/s7eA9F0lLBdKuL3/iawa5LeahqNzc3dxewkGKaWd5DJIVwMKzFccYxXSXHw38N3fiDWtbn0tJ9S1rT00rUJJJHZbi1QuViKFtgH72TkAE7uSeK6WigDyDw9+yj8OfDOpaDqFnYas97oE6zaTPda9fXDWCqCvkw75jshKsQYh8jDAYHAq9oH7M/w78M+LoPEWn6JNHeWs8t3aWsmo3MtjZzy582WC0aQwxO245KIOpxjJr1GigDy3w7+zP8PvC3ii01zT9Ku0ksbiS7sdPm1S6m0+xnfO+W3tHkMMTnc3KIMZOMZr1KiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigAooooAKKKKACiiigD/9k="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Прв пример - Поставување на парапет\n",
    "\n",
    "Првиот пример се однесува на проблем на поставување парапет (ограда) на правоаголна површина.\n",
    "\n",
    "Во внатрешноста на површината постои празен правоаголен простор со големина од $180m^{2}$. Задачата е да се најде најголемата можна оградена површина околу празниот простор, така што широчината на патеката за пешачење да е еднаква од сите страни, а да се искористат максимално $250m$ ограда.\n",
    "\n",
    "Графички проблемот е претставен на следната слика:\n",
    "\n",
    "![pr1.jpg](attachment:pr1.jpg)\n",
    "\n",
    "Од описот на проблемот следува дека како одлучувачки променливи може да се земат страните на правоаголната површина која што треба да се огради, односно:\n",
    "\n",
    "$X_{1}$ - должина на страната на правоаголната површина паралелна на кратката страна на празниот простор (10m), во единица мерка метри.\n",
    "<br>\n",
    "$X_{2}$ - должина на страната на правоаголната површина паралелна на долгата страна на празниот простор (18m), во единица мерка метри.\n",
    "\n",
    "Според главната цел на проблемот, што гласи дека треба да се најде најголемата површина која што може да се огради (без притоа да се нарушат ограничувањата), функцијата на цена на може да се претстави во следниот облик:\n",
    "\n",
    "$max f = X_{1}X_{2} - 180$, каде $X_{1}X_{2}$ е површината на оградениот правоаголник, додека со $180$ е претставена површината на внатрешниот празен простор.\n",
    "\n",
    "Ограничувањата на проблемот може да се претстават со следните равенки:\n",
    "\n",
    "$\\dfrac{(X_{1} - 10)}{2} = \\dfrac{(X_{2} - 18)}{2}$ - широчината на патеката за пешачење треба да биде еднаква од сите страни.\n",
    "\n",
    "$2\\cdot(X_{1} + X_{2}) + 2\\cdot(18 + 10) \\leq 250$ - најмногу може да се искористи $250m$ ограда.\n",
    "\n",
    "$X_{1}, X_{2} \\geq 0$ - одлучувачките променливи мора да бидат позитивни или еднакви на нула.\n",
    "\n",
    "По математичко средување на равенствата, се добива канонската форма на проблемот:\n",
    "\n",
    "<center>\n",
    "    <br>\n",
    "    $max f = X_{1}X_{2} - 180$\n",
    "    <br><br>\n",
    "    $0.5X_{1} - 0.5X_{2} + 4 \\leq 0$\n",
    "    <br>\n",
    "    $-0.5X_{1} + 0.5X_{2} - 4 \\leq 0$\n",
    "    <br>\n",
    "    $2X_{1} + 2X_{2} + 56 - 250 \\leq 0$\n",
    "</center>\n",
    "\n",
    "Со користење на овие равенства и со помош на функцијата *basinhopping* од *optimize* блокот на *SciPy* библиотеката, може со *Python* да се најде решението на овој проблем на оптимизација.\n",
    "\n",
    "*Basinhopping* алгоритамот се користи за наоѓање на глобалниот минимум на одредена функција на тој начин што комбинира алгоритам за глобално чекорење со алгоритам за одредување на локалниот минимум на секој чекор. Функцијата *minimize* од истата библиотека во некои одредени случаи ќе даде локален минимум како решение, а за разлика од неа *basinhopping* алгоритамот гарантира поголема сигурност за пронаоѓање на глобалниот минимум на функцијата која се оптимизира.\n",
    "\n",
    "Бидејќи алгоритамот кој што се користи врши минимизација на функција, потребно е проблемот да го запишеме во подруг облик, односно максимизирањето на функцијата на цена да премине во минимизирање, така што финалниот облик на проблемот е следниот:\n",
    "\n",
    "<center>\n",
    "    <br>\n",
    "    $min f = -X_{1}X_{2} + 180$\n",
    "    <br><br>\n",
    "    $0.5X_{1} - 0.5X_{2} + 4 \\leq 0$\n",
    "    <br>\n",
    "    $-0.5X_{1} + 0.5X_{2} - 4 \\leq 0$\n",
    "    <br>\n",
    "    $2X_{1} + 2X_{2} + 56 - 250 \\leq 0$\n",
    "</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.optimize import basinhopping, NonlinearConstraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_nlp(dcs_vars, cst_func, bnds, cons):\n",
    "    \n",
    "    # поставување на параметри на алгоритамот за минимизација\n",
    "    minimizer_args={'method': 'SLSQP', 'bounds': bnds, 'constraints': cons}\n",
    "    \n",
    "    # решавање на проблемот\n",
    "    solution = basinhopping(cst_func, dcs_vars, niter=100, minimizer_kwargs=minimizer_args, seed=157)\n",
    "\n",
    "    # добивање на пресметаните одлучувачки променливи\n",
    "    x = solution.x\n",
    "\n",
    "    # добивање на вредноста на функцијата на цена\n",
    "    cf_value = cst_func(x)\n",
    "    \n",
    "    return x, cf_value"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Оградената површина изнесува 2156.25 m², а ширината на патеката за пешачење е 17.25 m.\n",
      "\n",
      "Одлучувачки променливи:\n",
      "X₁ = 44.50 m\n",
      "X₂ = 52.50 m\n"
     ]
    }
   ],
   "source": [
    "# декларирање функција на цена\n",
    "def cost_function(x):\n",
    "    return -(x[0] * x[1] - 180)\n",
    "\n",
    "\n",
    "# декларирање на првото ограничување\n",
    "def constraint_1(x):\n",
    "    return (0.5 * x[0] - 0.5 * x[1] + 4)\n",
    "\n",
    "\n",
    "# декларирање на второто ограничување\n",
    "def constraint_2(x):\n",
    "    return (-0.5 * x[0] + 0.5 * x[1] - 4)\n",
    "\n",
    "\n",
    "# декларирање на третото ограничување\n",
    "def constraint_3(x):\n",
    "    return (2 * x[0] + 2 * x[1] + 56 - 250)\n",
    "\n",
    "\n",
    "# декларирање на одлучувачки променливи\n",
    "var = np.zeros(2)\n",
    "\n",
    "# поставување на ограничувања на променливите\n",
    "bounds = ((0, np.inf), (0, np.inf))\n",
    "\n",
    "# поставување на ограничувања на проблемот\n",
    "cons_1 = NonlinearConstraint(constraint_1, lb=-np.inf, ub=0.0)\n",
    "cons_2 = NonlinearConstraint(constraint_2, lb=-np.inf, ub=0.0)\n",
    "cons_3 = NonlinearConstraint(constraint_3, lb=-np.inf, ub=0.0)\n",
    "constraints = ([cons_1,cons_2, cons_3])\n",
    "\n",
    "decision_vars, cost_function_value = solve_nlp(dcs_vars=var, cst_func=cost_function, bnds=bounds, cons=constraints)\n",
    "\n",
    "# добивање на ширината на патеката за пешачење\n",
    "walkway_width = (decision_vars[0] - 10) / 2\n",
    "\n",
    "# принтање на вредноста на функцијата на цена\n",
    "print('Оградената површина изнесува {:.2f} m\\u00b2, а ширината на патеката за пешачење е {:.2f} m.\\n'.format(-cost_function_value, walkway_width))\n",
    "\n",
    "# принтање на вредноста на одлучувачките променливи\n",
    "print('Одлучувачки променливи:')\n",
    "print('X\\N{SUBSCRIPT ONE} = {:.2f} m'.format(decision_vars[0]))\n",
    "print('X\\N{SUBSCRIPT TWO} = {:.2f} m'.format(decision_vars[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Анализа на чувствителност\n",
    "\n",
    "Дополнително за конкретниот проблем, може да се направи и анализа на чувствителност.\n",
    "\n",
    "Единствената примена каде би имала смисла оваа постапка е во третото ограничување, односно до колку метри може да се намали дозволената должина на оградата:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Доколку се промени дозволената должина на ограда, основните променливи се следните: {x2: 52.5 - 0.25*Δ, x1: 44.5 - 0.25*Δ}\n",
      "\n",
      "Промената Δ која се појавува во ограничувањето мора да го задоволи следниот услов: Δ <= 178.0\n"
     ]
    }
   ],
   "source": [
    "import sympy as sym\n",
    "from sympy.solvers.inequalities import reduce_rational_inequalities\n",
    "\n",
    "d = sym.Symbol('\\u0394', real='True')\n",
    "x1 = sym.Symbol('x1', real='True')\n",
    "x2 = sym.Symbol('x2', real='True')\n",
    "\n",
    "solution_constraint_3 = sym.solve((0.5 * x1 - 0.5 * x2 + 4,\n",
    "                                  -0.5 * x1 + 0.5 * x2 - 4,\n",
    "                                   2 * x1 + 2 * x2 + 56 - 250 + d), (x1, x2, d))\n",
    "\n",
    "print('Доколку се промени дозволената должина на ограда, основните променливи се следните:', solution_constraint_3)\n",
    "\n",
    "constraint_3_range = reduce_rational_inequalities([[solution_constraint_3[x1] >= 0, solution_constraint_3[x2] >= 0]], d)\n",
    "\n",
    "print('\\nПромената \\u0394 која се појавува во ограничувањето мора да го задоволи следниот услов:', constraint_3_range)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Според добиените резултати, должината на оградата би можела да се намали до $250 - \\Delta = 72m$\n",
    "\n",
    "Доколку се промени ограничувањето според оваа величина, се добива следното решение на проблемот:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Оградената површина изнесува -180.00 m², а ширината на патеката за пешачење е -5.00 m.\n",
      "\n",
      "Одлучувачки променливи:\n",
      "X₁ = 0.00 m\n",
      "X₂ = 8.00 m\n"
     ]
    }
   ],
   "source": [
    "# декларирање функција на цена\n",
    "def cost_function(x):\n",
    "    return -(x[0] * x[1] - 180)\n",
    "\n",
    "\n",
    "# декларирање на првото ограничување\n",
    "def constraint_1(x):\n",
    "    return 0.5 * x[0] - 0.5 * x[1] + 4\n",
    "\n",
    "\n",
    "# декларирање на второто ограничување\n",
    "def constraint_2(x):\n",
    "    return -0.5 * x[0] + 0.5 * x[1] - 4\n",
    "\n",
    "\n",
    "# декларирање на третото ограничување\n",
    "def constraint_3(x):\n",
    "    return 2 * x[0] + 2 * x[1] + 56 - 72\n",
    "\n",
    "\n",
    "# декларирање на одлучувачки променливи\n",
    "var = np.zeros(2)\n",
    "\n",
    "# поставување на ограничувања на променливите\n",
    "bounds = ((0, np.inf), (0, np.inf))\n",
    "\n",
    "# поставување на ограничувања на проблемот\n",
    "cons_1 = NonlinearConstraint(constraint_1, lb=-np.inf, ub=0)\n",
    "cons_2 = NonlinearConstraint(constraint_2, lb=-np.inf, ub=0)\n",
    "cons_3 = NonlinearConstraint(constraint_3, lb=-np.inf, ub=0)\n",
    "constraints = ([cons_1,cons_2, cons_3])\n",
    "\n",
    "decision_vars, cost_function_value = solve_nlp(dcs_vars=var, cst_func=cost_function, bnds=bounds, cons=constraints)\n",
    "\n",
    "# добивање на ширината на патеката за пешачење\n",
    "walkway_width = (decision_vars[0] - 10) / 2\n",
    "\n",
    "# принтање на вредноста на функцијата на цена\n",
    "print('Оградената површина изнесува {:.2f} m\\u00b2, а ширината на патеката за пешачење е {:.2f} m.\\n'.format(-cost_function_value, walkway_width))\n",
    "\n",
    "# принтање на вредноста на одлучувачките променливи\n",
    "print('Одлучувачки променливи:')\n",
    "print('X\\N{SUBSCRIPT ONE} = {:.2f} m'.format(decision_vars[0]))\n",
    "print('X\\N{SUBSCRIPT TWO} = {:.2f} m'.format(decision_vars[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Може да се забележи дека, како што покажа и анализата на чувствителност, ова е минимално што може да се намали дозволената должина на оградата, бидејќи одлучувачката променлива $X_{1}$ добива вредност $0m$.\n",
    "\n",
    "Но, исто така може да се забележи дека добиената површина како и ширината на патеката за пешачење добива негативна вредност, што е невозможно. Ова е поради тоа што должините на страните на надворешниот правоаголник, односно одлучувачките променливи, се помали од должините на внатрешниот празен простор.\n",
    "\n",
    "За да се добие вистинската минимална дозволена должина на оградата, потребно е дополнително да се вметне и ограничувањето дека површината на надворешниот правоаголник мора да е поголема од површината на внатрешниот празен простор:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Доколку се промени дозволената должина на ограда, основните променливи се следните: [(-18.0000000000000, -10.0000000000000, 250.000000000000), (10.0000000000000, 18.0000000000000, 138.000000000000)]\n"
     ]
    }
   ],
   "source": [
    "import sympy as sym\n",
    "from sympy.solvers.inequalities import reduce_rational_inequalities\n",
    "\n",
    "d = sym.Symbol('\\u0394', real='True')\n",
    "x1 = sym.Symbol('x1', real='True')\n",
    "x2 = sym.Symbol('x2', real='True')\n",
    "\n",
    "solution_constraint_3 = sym.solve((0.5 * x1 - 0.5 * x2 + 4,\n",
    "                                  -0.5 * x1 + 0.5 * x2 - 4,\n",
    "                                   x1 * x2 - 180,\n",
    "                                   2 * x1 + 2 * x2 + 56 - 250 + d), (x1, x2, d))\n",
    "\n",
    "print('Доколку се промени дозволената должина на ограда, основните променливи се следните:', solution_constraint_3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Според добиените резултати постојат само две можни решенија, а тоа се:\n",
    "\n",
    "$X_{1} = -18m$\n",
    "<br>\n",
    "$X_{2} = -10m$\n",
    "<br>\n",
    "$\\Delta = 250$\n",
    "\n",
    "Ова решение е невозможно бидејќи $X_{1}, X_{2} \\geq 0$\n",
    "\n",
    "$X_{1} = 10m$\n",
    "<br>\n",
    "$X_{2} = 18m$\n",
    "<br>\n",
    "$\\Delta = 138$\n",
    "\n",
    "Според ова решение може да се воочи дека $X_{1}$ и $X_{2}$ се со иста должина како и страните на внатрешниот празен простор, што значи дека навистина е минималната дозволена должина на оградата, која би изнесувала $250 - \\Delta = 112m$, иако не е реалистично решение бидејќи во овој случај површината на патеката за пешачење би била $0m$. Следствено на ова, може да се заклучи дека $\\Delta \\leq 138$, односно најмалку може да се употреби $112m$ ограда."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Втор пример - Планирање на производство во фабрика за пијалоци\n",
    "\n",
    "Вториот пример во оваа проектна задача се однесува на проблем со планирање на производството во фабрика за пијалоци.\n",
    "\n",
    "Фабриката произведува два типа на пијалоци базирани на истиот сируп, *Coli* и *Pepsa*, на својата производствена линија *S1*.\n",
    "\n",
    "- По процесирањето, секој хектолитар на сируп произведува **40** единици со шишиња од 1.5l од пијалокот *Coli* и **20** единици со шишиња од 1.5l од пијалокот *Pepsa*. \n",
    "- Секој хектолитар на сируп ја кошта фирмата по **$150** за набавка и процесирање.\n",
    "- Производствената линија за пакување има капацитет за производство на **7,100** единици од 1.5l шишиња на секој час.\n",
    "- Фирмата работи **5 дена** во неделата по **8 часа** на ден.\n",
    "- Според неделната побарувачка, количеството на произведени единици од пијалокот *Pepsa* треба да биде барем половина од количеството на произведени единици од пијалокот *Coli*.\n",
    "\n",
    "Целта на решавање на овој проблем е да се добие оптималното дневно количество на произведени единици за секој од пијалоците, така што фирмата ќе го максимизира својот профит.\n",
    "\n",
    "Од описот на проблемот следува дека како одлучувачки променливи може да се земе бројот на произведени единици од поединечните пијалоци, односно:\n",
    "\n",
    "$X_{1}$ - број на произведени единици од пијалокот *Coli*.\n",
    "<br>\n",
    "$X_{2}$ - број на произведени единици од пијалокот *Pepsa*.\n",
    "\n",
    "Функцијата на цена е дадена во самиот проблем, односно таа претставува проценка на дневниот приход кој што го очекува фирмата минус потребните трошоци за набавување и процесирање на сирупот: \n",
    "\n",
    "\n",
    "<center>\n",
    "<br>\n",
    "$max f = 49000X_{1} - X_{1}^2 + 30X_{2} - 2X_{2}^2 - 150\\cdot(40X_{1} + 20X_{2})$\n",
    "</center>\n",
    "\n",
    "Ограничувањата на проблемот може да се претстават со следните равенки:\n",
    "\n",
    "$X_{1} + X_{2} \\leq 56800$ - вкупното количество на произведени единици мора да биде реално остварливо (7,100 единици на час, 8 часа на ден).\n",
    "\n",
    "$X_{2} \\geq \\dfrac{X_{1}}{2}$ - количеството на произведени единици од пијалокот *Pepsa* треба да биде барем половина од количеството на единици од пијалокот *Coli*.\n",
    "\n",
    "$X_{1}, X_{2} \\geq 0$ - одлучувачките променливи мора да бидат позитивни или еднакви на нула.\n",
    "\n",
    "По математичко средување на равенствата, се добива канонската форма на проблемот:\n",
    "\n",
    "<center>\n",
    "<br>\n",
    "$max f = 43000X_{1} - X_{1}^2 - 2970X_{2} - 2X_{2}^2$\n",
    "<br><br>\n",
    "$X_{1} + X_{2} - 56800 \\leq 0$\n",
    "<br>\n",
    "$X_{1} - 2X_{2} \\leq 0$\n",
    "</center>\n",
    "\n",
    "Бидејќи алгоритамот кој што се користи врши минимизација на функција, потребно е проблемот да го запишеме во подруг облик, односно максимизирањето на функцијата на цена да премине во минимизирање, така што финалниот облик на проблемот е следниот:\n",
    "\n",
    "<center>\n",
    "<br>\n",
    "$min f = -43000X_{1} + X_{1}^2 + 2970X_{2} + 2X_{2}^2$\n",
    "<br><br>\n",
    "$X_{1} + X_{2} - 56800 \\leq 0$\n",
    "<br>\n",
    "$X_{1} - 2X_{2} \\leq 0$\n",
    "</center>\n",
    "\n",
    "За решавање на проблемот повторно ја користиме функцијата *basinhopping*:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Остварениот профит на фирмата ќе изнесува: 287514496.15\n",
      "\n",
      "Одлучувачки променливи:\n",
      "X₁ = 13841.84 единици\n",
      "X₂ = 6912.26 единици\n"
     ]
    }
   ],
   "source": [
    "# декларирање функција на цена\n",
    "def cost_function(x):\n",
    "    return -43000*x[0] + x[0]**2 + 2970*x[1] + 2*x[1]**2\n",
    "\n",
    "\n",
    "# декларирање на првото ограничување\n",
    "def constraint_1(x):\n",
    "    return x[0] + x[1] - 56800\n",
    "\n",
    "\n",
    "# декларирање на второто ограничување\n",
    "def constraint_2(x):\n",
    "    return x[0] - 2 * x[1]\n",
    "\n",
    "\n",
    "# декларирање на одлучувачки променливи\n",
    "var = np.zeros(2)\n",
    "\n",
    "# поставување на ограничувања на променливите\n",
    "bounds = ((0, np.inf), (0, np.inf))\n",
    "\n",
    "# поставување на ограничувања на проблемот\n",
    "cons_1 = NonlinearConstraint(constraint_1, lb=-np.inf, ub=0)\n",
    "cons_2 = NonlinearConstraint(constraint_2, lb=-np.inf, ub=0)\n",
    "constraints = ([cons_1, cons_2])\n",
    "\n",
    "decision_vars, cost_function_value = solve_nlp(dcs_vars=var, cst_func=cost_function, bnds=bounds, cons=constraints)\n",
    "\n",
    "# принтање на вредноста на функцијата на цена\n",
    "print('Остварениот профит на фирмата ќе изнесува: {:.2f}\\n'.format(-cost_function_value))\n",
    "\n",
    "# принтање на вредноста на одлучувачките променливи\n",
    "print('Одлучувачки променливи:')\n",
    "print('X\\N{SUBSCRIPT ONE} = {:.2f} единици'.format(decision_vars[0]))\n",
    "print('X\\N{SUBSCRIPT TWO} = {:.2f} единици'.format(decision_vars[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Анализа на чувствителност\n",
    "\n",
    "За овој проблем, анализата на чувствителност ќе се изведе со користење на ***Kuhn-Tucker*** условите. \n",
    "\n",
    "Преку оваа постапка може целосно да се реши проблемот, без користење на апроксимирачки алгоритми што се користеа до сега. Покрај тоа, на овој начин се добиваат и дополнителни информации за проблемот кои што се корисни за анализата на чувствителност.\n",
    "\n",
    "Постапката е базирана на користење на Лагранжова функција и поставување на т.н. *Kuhn-Tucker*-ови услови кои мораат да бидат исполнети. На крај решавањето се сведува на решавање на систем равенки, што не е тешко да се изведе."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import sympy as sym\n",
    "\n",
    "l1 = sym.Symbol('\\u03BB_1', real='True')\n",
    "l2 = sym.Symbol('\\u03BB_2', real='True')\n",
    "x1 = sym.Symbol('x1', real='True')\n",
    "x2 = sym.Symbol('x2', real='True')\n",
    "\n",
    "# декларирање на Лагранжовата функција, што претставува збир од функцијата на цена со секој од условите помножени \n",
    "# со соодветен Лагранжов мултипликатор\n",
    "L = -43000*x1 + x1**2 + 2970*x2 + 2*x2**2 + l1*(x1 + x2 - 56800) + l2*(x1 - 2 * x2)\n",
    "\n",
    "# услови на Kuhn-Tucker\n",
    "\n",
    "# услов на градиент\n",
    "dLx1 = L.diff(x1)\n",
    "dLx2 = L.diff(x2)\n",
    "\n",
    "# услов на изводливост\n",
    "f1 = x1 + x2 - 56800\n",
    "f2 = x1 - 2*x2\n",
    "\n",
    "# услов на ортогоналност\n",
    "o1 = l1*(x1 + x2 - 56800)\n",
    "o2 = l2*(x1 - 2*x2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "За да се реши проблемот, треба да се претпостават вредностите на Лагранжовите мултипликатори:\n",
    "\n",
    "**1)** $\\lambda_1 = \\lambda_2 = 0$\n",
    "\n",
    "Од условот на градиент следува:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x1 = 21500.0\n",
      "x2 = -742.5\n"
     ]
    }
   ],
   "source": [
    "var_1 = sym.solve(dLx1.subs({l1: 0, l2: 0}), x1)  \n",
    "var_2 = sym.solve(dLx2.subs({l1: 0, l2: 0}), x2)\n",
    "\n",
    "print('x1 = {}\\nx2 = {}'.format(float(var_1[0]), float(var_2[0])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Од овие вредности може да се забележи дека променливата $X_{1}$ е во ред, но променливата $X_{2}$ добива негативна вредност, што не е дозволено.\n",
    "\n",
    "Поради ова, наредно ќе се претпостават следните вредности за Лагранжовите мултипликатори:\n",
    "\n",
    "**2)** $\\lambda_1 = 0; \\lambda_2 \\neq 0$\n",
    "\n",
    "Од условот на градиент следува:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 109,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x1 = 21500 - λ_2/2\n",
      "x2 = λ_2/2 - 1485/2\n"
     ]
    }
   ],
   "source": [
    "eq_1 = sym.solve(dLx1.subs(l1, 0), x1)[0]\n",
    "eq_2 = sym.solve(dLx2.subs(l1, 0), x2)[0]\n",
    "\n",
    "print('x1 = {}\\nx2 = {}'.format(eq_1, eq_2))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Може да се воочи дека сега како резултат на замената се добиваат две равенки со три непознати. \n",
    "\n",
    "За да може да се најдат решенијата на равенките, потребно е да се добие уште една равенка која ги содржи двете променливи. Тоа може да се оствари со заменување на претпоставените вредности за Лагранжовите мултипликатори во останатите услови на *Kuhn-Tucker*.\n",
    "\n",
    "Бидејќи $\\lambda_{1} = 0$, единствена равенка која што може да се искористи е вториот ортогонален услов, односно:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(41515/3, 41515/6, 45970/3), (21500, -1485/2, 0)]\n"
     ]
    }
   ],
   "source": [
    "eq_3 = o2\n",
    "\n",
    "solution = sym.solve((eq_1 - x1, eq_2 - x2, eq_3 - 0), (x1, x2, l2))\n",
    "\n",
    "print(solution)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Бидејќи во едно од можните решенија променливата $X_{2}$ има негативна вредност, решението се зема за невозможно, па вистинското решение на проблемот е:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X₁ = 13838.333333333334\n",
      "X₂ = 6919.166666666667\n",
      "λ₂ = 15323.333333333334\n"
     ]
    }
   ],
   "source": [
    "x1_s, x2_s, l2_s = solution[0]\n",
    "\n",
    "print('X\\N{SUBSCRIPT ONE} = {}\\nX\\N{SUBSCRIPT TWO} = {}\\n\\u03BB\\N{SUBSCRIPT TWO} = {}'.format(float(x1_s), float(x2_s), float(l2_s)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Како што може да се забележи, резултатите се скоро еднакви со оние добиени со апроксимирачко минимизирање на функцијата на цена користејќи го *basinhopping* алгоритамот.\n",
    "\n",
    "За разлика од претходниот метод, овде имаме и податок за Лагранжовите мултипликатори $\\lambda_{1}$ и $\\lambda_{2}$, што не претставуваат ништо друго туку скриените цени или *shadow prices* на ограничувањата во проблемот.\n",
    "\n",
    "Според тоа, може да заклучиме дека бидејќи $\\lambda_{1} = 0$, фирмата нема да има корист од додавање на повеќе работни часови во денот, бидејќи не се искористени ни тие што веќе ги има на располагање. Тоа може да се потврди и со едноставна проверка:\n",
    "\n",
    "<center>\n",
    "    <br>\n",
    "    $X_{1} + X_{2} = 13838 + 6919 = 20757 < 56800$\n",
    "</center>\n",
    "\n",
    "Од $\\lambda_{2} = 15323$ може да се увиди дека за секоја дополнителна произведена единица од пијалокот *Coli*, дневниот профит би се зголемил за ***$15,323***."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}