MILP_with_Solver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1be7e4ca",
   "metadata": {},
   "source": [
    "## We start be loading some modules and helpers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32deec5f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "sys.path.append('.')\n",
    "import numpy as np\n",
    "from scipy import linalg\n",
    "from helpers import *\n",
    "\n",
    "from smt.sampling_methods import LHS\n",
    "\n",
    "import scipy.interpolate\n",
    "\n",
    "from docplex.mp.model import Model\n",
    "\n",
    "from qiskit_optimization.algorithms import CplexOptimizer,CobylaOptimizer\n",
    "from qiskit_optimization.translators import from_docplex_mp\n",
    "\n",
    "from qiskit.circuit.library import EfficientSU2\n",
    "from qiskit import Aer\n",
    "from qiskit.algorithms import VQE\n",
    "from qiskit_optimization.algorithms import MinimumEigenOptimizer\n",
    "from qiskit.algorithms.optimizers import COBYLA, L_BFGS_B, GSLS, GradientDescent\n",
    "\n",
    "from qiskit.circuit.library import RealAmplitudes\n",
    "from qiskit import QuantumCircuit, ClassicalRegister, QuantumRegister\n",
    "from qiskit.circuit import ParameterVector\n",
    "from scipy import optimize\n",
    "\n",
    "from scipy.optimize import minimize\n",
    "from smt.surrogate_models import KRG\n",
    "\n",
    "try:\n",
    "    import matplotlib.pyplot as plt\n",
    "    plot_status = True\n",
    "except:\n",
    "    plot_status = False\n",
    "from matplotlib import cm"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b66cead",
   "metadata": {},
   "source": [
    "##  Primary classes defined in the following cells are:\n",
    "- BinaryProblemBody: represents a pure QUBO problem via (objective function includes constraints via penalty term as soft constraints whose continuous variables have to be fixed for evaluation on the quantum device)\n",
    "- VQE_Optimizer: QUBO optimizer using VQE. The ansatz function can be choose by preference (we offer three different ansatz functions)\n",
    "- MILPproblem: a class representing the entire MILP (or here a MBLP) problem including objective function, constraints and penalty\n",
    "### Solvers\n",
    "We implemented two different strategies to solve the problem. \n",
    "- MILPAlternatingsolver: solve the MILPproblem with alternating approach where a QUBO with fixed continuous variables and a continuous problem with fixed binary variables are optimized in an alternating fashion.\n",
    "- MILPSurrogateSovler: solve the MILPproblem by building a purely continuous surrogate model and optimizing and refining the surrogate model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "173cfea4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# BinaryProblemBody requires two initial parameter: the objective function Expr (callable) and\n",
    "# nBinary, the number of binary variables\n",
    "# It uses the qiskit min_eigen_optimizer to build a QUBO problem\n",
    "class BinaryProblemBody:\n",
    "    def __init__(self,Expr,nBinary):\n",
    "        self.md = Model(\"MixedBinaryContinous\")\n",
    "        self.x = []\n",
    "        for i in range(nBinary):\n",
    "            self.x.append(self.md.binary_var(\"x\"+str(i)))\n",
    "        self.initExpr = Expr\n",
    "        self.opti = []\n",
    "        self.optivalue = 1000\n",
    "\n",
    "    def updateBodyExpr(self,EQConstrain):\n",
    "        self.Expr = lambda u,x: self.initExpr(u,x) + sum([con(u,x) for con in EQConstrain])\n",
    "\n",
    "    def cleanConstrain(self):\n",
    "        self.md.clear_constraints()\n",
    "        \n",
    "    def UpdateExpr(self,u):\n",
    "        # evaluate the expression with fixed continous value and\n",
    "        # update the minimization problem for the binary model. \n",
    "        self.md.minimize(self.Expr(u,self.x))\n",
    "        \n",
    "    def EvaNewBinaryOpti(self,min_eigen_optimizer):        \n",
    "        # Generate QUBO from Model md and solve it\n",
    "        qp = from_docplex_mp(self.md)\n",
    "        self.opti = min_eigen_optimizer.solve(qp).x\n",
    "        self.optivalue = min_eigen_optimizer.solve(qp).fval"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "55abab73",
   "metadata": {},
   "outputs": [],
   "source": [
    "# VQE_Optimizer requires only number of quibits to instantiate the class.\n",
    "class VQE_Optimizer:\n",
    "    def __init__(self,nqubits,WaveFunType = \"RealAmplitudes\", backend=Aer.get_backend('aer_simulator_statevector'), \n",
    "                 optimizer=COBYLA()):\n",
    "        self.backend = backend\n",
    "        self.nqubits = nqubits\n",
    "        self.optimizer = optimizer\n",
    "        self.circuit_plot = None\n",
    "        self.__generateOptimizer(WaveFunType)\n",
    "    \n",
    "    # Define an VQE optimizer, i.e. define an ansatz and an instance of qiskit's VQE algorithm\n",
    "    def __generateOptimizer(self,WaveFunType):\n",
    "        ## part 1: manually set the wave up\n",
    "        if(WaveFunType==\"manual\"):\n",
    "            params =ParameterVector('theta', length=self.nqubits*2)\n",
    "            it = iter(params)\n",
    "            wavefunction = QuantumCircuit(self.nqubits)\n",
    "            for i in range(0, self.nqubits):\n",
    "                wavefunction.ry(2 * next(it), i)\n",
    "\n",
    "            #adding extra layers\n",
    "            for i in range(0,self.nqubits-1):\n",
    "                wavefunction.cz(0,1)\n",
    "\n",
    "                for j in range(0, self.nqubits):\n",
    "                    wavefunction.ry(2 * next(it), j)\n",
    "\n",
    "        ## part 2:   EfficientSU2 ansatz\n",
    "        elif(WaveFunType ==\"EfficientSU2\"):\n",
    "            wavefunction = EfficientSU2(num_qubits = self.nqubits, entanglement=\"pairwise\")\n",
    "\n",
    "        ## part 3:   RealAmplitudes ansatz\n",
    "        elif(WaveFunType ==\"RealAmplitudes\"):\n",
    "            wavefunction = RealAmplitudes(num_qubits=self.nqubits, reps=3)\n",
    "        else:\n",
    "            print(\"the given ansatz type is currently not supported\")\n",
    "        print(wavefunction)\n",
    "        self.circuit_plot = wavefunction.decompose().draw(output='mpl')\n",
    "      \n",
    "        #######using VQE######\n",
    "        vqe = VQE(ansatz=wavefunction, optimizer=self.optimizer, quantum_instance=self.backend)\n",
    "        self.min_eigen_optimizer = MinimumEigenOptimizer(vqe)\n",
    "        ######################"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70029d83",
   "metadata": {},
   "outputs": [],
   "source": [
    "# MILPproblem requires BinaryProblemBody, VQE_optimizer and Penalty value \n",
    "# as well as number of classical and binary variables and boundary for continuous variables\n",
    "\n",
    "# The add_constrain method can add constrains from outside of the class and the BinaryProblemBody's \n",
    "# objective function via the updateConstrain method by adding a penalty term. \n",
    "# Currently, it supports for equality constrains only.\n",
    "\n",
    "class MILPproblem:\n",
    "    def __init__(self,body,VQE_Optimizer,Penalty, num_cl_vars, num_qbits, bounds):\n",
    "        self.num_cl_vars = num_cl_vars\n",
    "        self.num_qbits = num_qbits\n",
    "        self.bounds = bounds\n",
    "        self.body = body\n",
    "        self.Optimizer = VQE_Optimizer\n",
    "        self.Constrains = []\n",
    "        self.Penalty = Penalty\n",
    "        \n",
    "    def CalcBinaryOptimum(self):\n",
    "        self.body.EvaNewBinaryOpti(self.Optimizer.min_eigen_optimizer)\n",
    "        \n",
    "    def continousProblem(self,u,x):\n",
    "        #Two parts are involved, but x(binary parts) will be treated as parameter(fixed) \n",
    "        return self.body.Expr(u,x)\n",
    "        \n",
    "    def updateBinary(self,u):\n",
    "        self.body.UpdateExpr(u)\n",
    "        \n",
    "    def getBinaryOpti(self):\n",
    "        return self.body.opti\n",
    "    \n",
    "    def getBinaryOptiValue(self):\n",
    "        return self.body.optivalue\n",
    "    \n",
    "    def add_constrain(self,Expr,Op,target):\n",
    "        self.Constrains.append(constrains(Expr,Op,target))\n",
    "    \n",
    "    def updateConstrain(self):\n",
    "        self.EQConstrain = []\n",
    "        for cons in self.Constrains:\n",
    "            if(cons.op == \"==\"):\n",
    "                #update the subject function(Expr) according to the equality constrain\n",
    "                # Expr = Expr + penalty * (constrain-target)^2\n",
    "                self.EQConstrain.append(lambda u,x: self.Penalty*pow((cons.Expr(u,x)-cons.Target),2))\n",
    "        self.body.updateBodyExpr(self.EQConstrain)\n",
    "        \n",
    "    def removeConstrain(self):\n",
    "        self.body.cleanConstrain()\n",
    "        \n",
    "    def testVadility(self,res):\n",
    "        delta = 1e-04\n",
    "        for cons in self.EQConstrain:\n",
    "            if(cons(res[0],res[1])>delta):\n",
    "                return False\n",
    "        return True\n",
    "    \n",
    "    def getBinProbMinObjValue(self,u):\n",
    "        self.body.UpdateExpr(u)\n",
    "        self.CalcBinaryOptimum()\n",
    "        return self.body.optivalue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "afd1ea5d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# The MILPAlternatingsolver is initialized with a MILPproblem and solves it with an alternating strategy.\n",
    "# The alternating strategy works as following: \n",
    "#    1. Give an initial guess of the continous variables\n",
    "#    2. Update the continous part of the body expression of the MILPproblem, so that it will be a pure QUBO problem; solve the\n",
    "#       QUBO problem via the VQE solver.\n",
    "#    3. Update the body expression of the MILPproblem via the result of step 2. \n",
    "#       The problem becomes pure continous, it will be solved with an classical optimizer.\n",
    "#    4. If the generated sequence starts 'converging' (i.e. update between old and new small), then exit, else go to step 2.\n",
    "class MILPAlternatingsolver:\n",
    "    def __init__(self, MILPbody):\n",
    "        self.MILPbody = MILPbody\n",
    "        \n",
    "    def __distance(self,x,x1):\n",
    "        x_ = np.array(x)\n",
    "        x1_ = np.array(x1)\n",
    "        return np.linalg.norm(x_-x1_,2)\n",
    "    \n",
    "    def solver(self,u_init, verbose=False):\n",
    "        bounds = self.MILPbody.bounds\n",
    "        solution_old = [0]*(self.MILPbody.num_cl_vars+self.MILPbody.num_qbits)\n",
    "        step = 0\n",
    "        self.MILPbody.updateBinary(u_init)\n",
    "        # use alternating strategy until update between old/new value becomes small (or until at most 1e4 steps done)\n",
    "        while step<1E4:\n",
    "            step = step + 1\n",
    "            #find optimum of binary problem\n",
    "            self.MILPbody.CalcBinaryOptimum()\n",
    "            x =self.MILPbody.getBinaryOpti()\n",
    "\n",
    "            # Get constrains with fixed binary value\n",
    "            constr = []\n",
    "            for constraint in self.MILPbody.EQConstrain:\n",
    "                constr.append({'type': 'eq', 'fun': lambda u: constraint(u,x)})\n",
    "            # Optimize purely continuous problem\n",
    "            res = optimize.shgo(self.MILPbody.body.initExpr, bounds, args = (x,), constraints=constr)\n",
    "            u = res.x\n",
    "            # update binary problem\n",
    "            self.MILPbody.updateBinary(u)\n",
    "\n",
    "            solution_new = list(x) + list(u)\n",
    "            if verbose:\n",
    "                displayDetail(step,u,x)\n",
    "            if(self.__distance(solution_new,solution_old)<1E-8):\n",
    "                break\n",
    "            solution_old = solution_new\n",
    "        return [u,x] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a54b2ed5",
   "metadata": {},
   "outputs": [],
   "source": [
    "class MILPSurrogateSolver:\n",
    "    def __init__(self, MILPbody):\n",
    "        self.MILPbody = MILPbody\n",
    "        \n",
    "    def __SBO(self,point):\n",
    "        return self.t.predict_values(point)\n",
    "    \n",
    "    def __findCurrentOptimum(self,obj,npoints,bound):\n",
    "        r_TR_or = [bound[i][1]-bound[i][0] for i in range(len(bound))]\n",
    "        # generate initial points for optimization\n",
    "        x_start = np.zeros([80,len(bound)])\n",
    "        for i in range(len(bound)):\n",
    "            x_start[:,i] = bound[i][0] + np.random.rand(80)*r_TR_or[i]\n",
    "        #define the objective function:\n",
    "        opt_all = np.array([minimize(lambda x: float(obj(x)), x_st, method='SLSQP', bounds=bound) for x_st in x_start])\n",
    "        opt_success = opt_all[[opt_i['success'] for opt_i in opt_all]]\n",
    "        obj_success = np.array([opt_i['fun'] for opt_i in opt_success])\n",
    "        #return the optimum point\n",
    "        ind_min = np.argmin(obj_success)\n",
    "        opt = opt_success[ind_min]\n",
    "        return opt['x']\n",
    "    \n",
    "    def solver(self,npoints,bound,maxIter):\n",
    "        # generate initial sampling points and calculate their corresponding minmum value\n",
    "        sampling = LHS(xlimits = bound, criterion='ese', random_state=1)\n",
    "        # u_base and y_base are used as the data set for training the surrogate model\n",
    "        self.u_base = sampling(npoints)\n",
    "        self.y_base = []\n",
    "        # y_old and x_old are used to record the best point of surrogate model for last step\n",
    "        x_old  = np.array([-1000,-1000])\n",
    "        y_old  = 1\n",
    "        \n",
    "        # evaluate all initial data points in base data set in order to build the surrogate model\n",
    "        for u in self.u_base:\n",
    "            self.y_base.append(self.MILPbody.getBinProbMinObjValue(u))\n",
    "        # building surrogate using Kriging model\n",
    "        # more information please refer:\n",
    "        # https://www.sciencedirect.com/science/article/pii/S100093611930041X\n",
    "        self.t = KRG(theta0=[1e-4]*len(bound),print_prediction = False)\n",
    "        self.t.set_training_values(self.u_base,np.array(self.y_base))\n",
    "        self.t.train()\n",
    "        \n",
    "        # this objective funtion stands for the finding the best point of the current surrogate model\n",
    "        # which is a pure continuous problem. \n",
    "        obj_k = lambda x: self.__SBO(np.atleast_2d(x))\n",
    "        \n",
    "        for i in range(maxIter):# start to refine the surrogate model by adding maximal 'maxInter' points \n",
    "            #find the best point of the current surrogate model\n",
    "            x_et_k = self.__findCurrentOptimum(obj_k,80,bound)\n",
    "            \n",
    "            #check if the result is converged, if yes, return the result.\n",
    "            diff = x_et_k-x_old\n",
    "            if(np.linalg.norm(diff)<1e-3):\n",
    "                print(\"find optimal and break\")\n",
    "                self.MILPbody.updateBinary(x_et_k)\n",
    "                self.MILPbody.CalcBinaryOptimum()\n",
    "                x_binary =self.MILPbody.getBinaryOpti()\n",
    "                return [list(x_et_k),list(x_binary)]\n",
    "            \n",
    "            #evaluate the new data\n",
    "            #y_et_k is the true value calculated based on the current best point position.\n",
    "            y_et_k = self.MILPbody.getBinProbMinObjValue(x_et_k)\n",
    "            if(y_et_k<y_old):\n",
    "                center = x_et_k\n",
    "            else:\n",
    "                center = x_old\n",
    "            #converting nparray into list\n",
    "            self.y_base.append(y_et_k)\n",
    "            ulist = self.u_base.tolist()\n",
    "            ulist.append(x_et_k.tolist())\n",
    "            self.u_base = np.array(ulist)\n",
    "            x_old = x_et_k\n",
    "            y_old = y_et_k\n",
    "            \n",
    "            # adding curent best point to the base data site\n",
    "            self.t.set_training_values(self.u_base,np.array(self.y_base))\n",
    "            self.t.train()\n",
    "        self.MILPbody.updateBinary(center)\n",
    "        self.MILPbody.CalcBinaryOptimum()\n",
    "        x_binary =self.MILPbody.getBinaryOpti()\n",
    "        return [list(center),list(x_binary)]\n",
    "    \n",
    "    def visualization(self,bound):\n",
    "        xlabel = np.linspace(bound[0][0], bound[0][1], 50)\n",
    "        ylabel = np.linspace(bound[1][0], bound[1][1],50)\n",
    "        resSM = []\n",
    "        varSM = []\n",
    "        x2D = []\n",
    "        for xx0 in xlabel:\n",
    "            for xx1 in ylabel:\n",
    "                x2D.append(np.array([xx0,xx1]))\n",
    "                resSM.append(self.t.predict_values(np.array([[xx0,xx1]])))\n",
    "                varSM.append(self.t.predict_variances(np.array([[xx0,xx1]])))\n",
    "        resSM = np.array(resSM)\n",
    "        resSM = resSM.reshape((50,50)).T\n",
    "        varSM = np.array(varSM)\n",
    "        varSM = varSM.reshape((50,50)).T                     \n",
    "        X,Y = np.meshgrid(xlabel,ylabel)\n",
    "\n",
    "        fig = plt.figure(figsize=(15, 10))\n",
    "        ax = fig.gca(projection='3d')\n",
    "        ax.scatter(self.u_base[:,0], self.u_base[:,1], self.y_base, zdir='z', marker='x', c='b', s=200, label='DOE')\n",
    "        surf = ax.plot_surface(X, Y, resSM, cmap=cm.coolwarm,\n",
    "                               linewidth=0, antialiased=False,alpha=0.5)\n",
    "\n",
    "        plt.legend()\n",
    "        plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6dcb9ce7",
   "metadata": {},
   "source": [
    "# Creating a MILP problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "dc32fe6d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# We define the following MILP (binary)\n",
    "#     objektive function: x1*u1 - x1*x2 + x1*x2*u2 - u1+ 5*u2 + u1*u2 + u1*u1 + u2*u2\n",
    "#     binary variables : x1,x2\n",
    "#     continous variables: u1,u2\n",
    "#     constraint (bounds) for continuous variables: -2<u1<2,-2<u2<3\n",
    "#     and some additional constraints, see below.\n",
    "\n",
    "# We have to specify the number of qubits (i.e. binary variables)\n",
    "nQubit = 2\n",
    "# and the number of continuous variables\n",
    "num_classical_vars = 2\n",
    "\n",
    "# the objective function\n",
    "BodyExpr = lambda u,x:x[0]*u[0] - x[1]*u[1] + x[0]*x[1]\\\n",
    "                    - u[0]+ 5*u[1]+ u[0]*u[1]+u[0]*u[0]\n",
    "\n",
    "# and we define the additional constraints\n",
    "Constrain1 = lambda u,x: x[0]*u[0] - x[1]*u[1]\n",
    "Constrain2 = lambda u,x: x[0]*u[0] + x[1]*u[1]\n",
    "# those will be implemented as soft-constraints in the QUBO with penalty factor, i.e. something like\n",
    "# ObjectiveFunction = lambda u,x:BodyExpr(u,x)+Penalty*pow(Constrain1(u,x),2)\n",
    "Penalty = 30\n",
    "# and bounds\n",
    "bounds_classical_vars = np.array([(-2,2),(-2,3)])\n",
    "\n",
    "# create binary optimization problem with the expression above\n",
    "binaryP = BinaryProblemBody(BodyExpr,nQubit)\n",
    "\n",
    "# then we set up the VQE algorithm with ansatz \"RealAmplitudes\" \n",
    "backend = Aer.get_backend('qasm_simulator') # We explicitly choose a backend; otherwise default aer_simulator_statevector will be used\n",
    "optimizer = COBYLA() # default if not set; other options: L_BFGS_B, GSLS, GradientDescent\n",
    "BinaryOptimizer = VQE_Optimizer(nQubit,\"RealAmplitudes\", backend, optimizer=optimizer)   \n",
    "\n",
    "# and finally we create the MILP problem\n",
    "MILP = MILPproblem(binaryP,BinaryOptimizer,Penalty, nQubit, num_classical_vars, bounds_classical_vars)\n",
    "\n",
    "# and add one constraint MILPproblem in the sense Constrain1 = 2.0 (last argument specifies target value for the constraint)\n",
    "MILP.add_constrain(Constrain1,\"==\",2.0)\n",
    "# MILP.add_constrain(Constrain2,\"==\",2.0)\n",
    "\n",
    "# after adding, the constrain should be updated. The \n",
    "MILP.updateConstrain()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a13225b",
   "metadata": {},
   "source": [
    "# Using the alternating solver to find the solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a733628b",
   "metadata": {},
   "outputs": [],
   "source": [
    "solver = MILPAlternatingsolver(MILP)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "899ebbea",
   "metadata": {},
   "source": [
    "## The final optimum is initial condition depedent\n",
    "An expected problem is that due to local minima or flat areas we might not find a global optimum. See the following example using onle one constraint and to different initial values for the continuous variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "99439188",
   "metadata": {},
   "outputs": [],
   "source": [
    "res1 = solver.solver([1,0], verbose=True)\n",
    "pretty_print_optimizer(res1)\n",
    "print(f\"Optimum: {BodyExpr(res1[0],res1[1])}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b714f48e",
   "metadata": {},
   "source": [
    "For another inital value we arrive at another point (and we reduce the solvers verbosity)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5afc7292",
   "metadata": {},
   "outputs": [],
   "source": [
    "res2 = solver.solver([0,-2])\n",
    "pretty_print_optimizer(res2)\n",
    "print(f\"Optimum: {BodyExpr(res2[0],res2[1])}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "29fa1e63",
   "metadata": {},
   "source": [
    "## Use LHS to find the possibly best optimum\n",
    "One way to deal with this problem is to use many different intial values for the continuous variables and take the minimum of the resulting local 'optima'. We use latin hypercube sampling (LHS) to define the initial values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "664e8feb",
   "metadata": {},
   "outputs": [],
   "source": [
    "#using LHS to generate the initial points\n",
    "#select the best result starting from different initial points\n",
    "#it is not guaranteed to find the global optimum, but works well for a lot situation\n",
    "node = 10\n",
    "sampling = LHS(xlimits = bounds_classical_vars, criterion='ese', random_state=1)\n",
    "Uinit = sampling(node)\n",
    "optimum =float('Inf')\n",
    "for uinit in Uinit:\n",
    "    resi = solver.solver(uinit)\n",
    "    #check if the solution is valid (constraints satisfied) and smaller than the old optimum\n",
    "    if(MILP.testVadility(resi) and BodyExpr(resi[0],resi[1])<optimum):\n",
    "        optimum = BodyExpr(resi[0],resi[1])\n",
    "        optimumResult = resi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eca2bfcd",
   "metadata": {},
   "outputs": [],
   "source": [
    "if optimum<float('Inf'):\n",
    "    print(f\"Optimum found is: {optimum} with\")\n",
    "    pretty_print_optimizer(optimumResult)\n",
    "else:\n",
    "    print(\"No feasible optimum found.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "274b6a43",
   "metadata": {},
   "source": [
    "# Using the surrogate model to solve the problem\n",
    "As an alternative approach we can use the surrogate solver that builds and refines a continuous Kriging surrogate model by sampling the MILP problem, i.e. optimizing the binary problem for fixed continuous sampling points.\n",
    "Then the surrogate model is minimized and the resulting point is added as sampling point if the update is larger than 1e-3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f3c20cc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "SurrogateSolver = MILPSurrogateSolver(MILP)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "08907bda",
   "metadata": {},
   "outputs": [],
   "source": [
    "# the solver of surrogate mode based method requires 3 input paramters:\n",
    "# 1. Number of initial point of the surrogate model\n",
    "# 2. Boundary of the continuous variable\n",
    "# 3. Maximum iteration of refinement\n",
    "center = SurrogateSolver.solver(40,bounds_classical_vars,20)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "47b960fe",
   "metadata": {},
   "source": [
    "We can test if the solution is a feasible solution (since the surrogate model does not necessarily satisfies the constraints due to the formulation as soft constraints we should check):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0aa485a",
   "metadata": {},
   "outputs": [],
   "source": [
    "MILP.testVadility(center)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8f3078d",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Optimum found: {BodyExpr(center[0],center[1])}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8cc847fd",
   "metadata": {},
   "source": [
    "Let's visualize how the surrogate model looked like and how the refinement worked"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "720de027",
   "metadata": {},
   "outputs": [],
   "source": [
    "SurrogateSolver.visualization(bounds_classical_vars)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "529b87dc",
   "metadata": {},
   "source": [
    "## Let's add another constraint"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "efcffd7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Recall that we considered the problem\n",
    "#     objective function: x1*u1 - x1*x2 + x1*x2*u2 - u1+ 5*u2 + u1*u2 + u1*u1 + u2*u2\n",
    "#     binary variables: x1,x2\n",
    "#     continous variables: u1,u2\n",
    "#     with constraints: -2<u1<2,\n",
    "#                       -2<u2<3,\n",
    "#                       x[0]*u[0] - x[1]*u[1] = 2.0\n",
    "\n",
    "# We add the constraint \n",
    "# Constrain2: x[0]*u[0] + x[1]*u[1] = 2.0\n",
    "MILP.add_constrain(Constrain2,\"==\",2.0)\n",
    "\n",
    "# after adding, the constrain should be updated.\n",
    "MILP.updateConstrain()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f0418a8",
   "metadata": {},
   "outputs": [],
   "source": [
    "node = 10\n",
    "sampling = LHS(xlimits = bounds_classical_vars, criterion='ese', random_state=1)\n",
    "Uinit = sampling(node)\n",
    "solver = MILPAlternatingsolver(MILP)\n",
    "optimum =float(\"Inf\")\n",
    "for uinit in Uinit:\n",
    "    resi = solver.solver(uinit)\n",
    "    #check if the solution is valid (constraints satisfied) and smaller than the old optimum\n",
    "    if(MILP.testVadility(resi) and BodyExpr(resi[0],resi[1])<optimum):\n",
    "        optimum = BodyExpr(resi[0],resi[1])\n",
    "        optimumResult = resi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c20a74b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "if optimum<float('Inf'):\n",
    "    print(f\"Optimum found is: {optimum} with\")\n",
    "    pretty_print_optimizer(optimumResult)\n",
    "else:\n",
    "    print(\"No feasible optimum found.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6ee16fb6",
   "metadata": {},
   "outputs": [],
   "source": [
    "SurrogateSolver = MILPSurrogateSolver(MILP)\n",
    "center = SurrogateSolver.solver(40,bounds_classical_vars,20)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e3bb27fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "pretty_print_optimizer(center)\n",
    "print(BodyExpr(center[0],center[1]))\n",
    "print(MILP.testVadility(center))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef4ab717",
   "metadata": {},
   "source": [
    "In the first example, alternating solver showed a good balance between accuracy and calculation cost. However, in the second \n",
    "example, it is failed to find a feasible point."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "eaa7e234",
   "metadata": {},
   "source": [
    "## Conclusions\n",
    "Both algorithms delivered in the test cases feasible results. Since we needed less (expensive) QUBO evaluations for the alternating directions approach might be better suited for current small NISQ devices. However, in particular for larger problems and on future quantum devices, the surrogate model might also be interesting since it is not purely local as the alternating directions approach."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "id": "382afcf6",
   "metadata": {},
   "source": [
    "**Copyright 2023 Tobias Haas and Qifeng Pan, University of Stuttgart - HLRS**\n",
    "\n",
    "Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "you may not use this file except in compliance with the License.\n",
    "You may obtain a copy of the License at\n",
    "\n",
    "    http://www.apache.org/licenses/LICENSE-2.0\n",
    "\n",
    "Unless required by applicable law or agreed to in writing, software\n",
    "distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "See the License for the specific language governing permissions and\n",
    "limitations under the License."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.9.13 (main, Aug 25 2022, 23:51:50) [MSC v.1916 64 bit (AMD64)]"
  },
  "vscode": {
   "interpreter": {
    "hash": "22af032344c768e1d157e1b38d6a78b920a14047785a8ab4488015d1714ac992"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}