Skip to content
This repository has been archived by the owner on Dec 4, 2024. It is now read-only.

Commit

Permalink
Update 1019
Browse files Browse the repository at this point in the history
  • Loading branch information
@Dran-Z
Dran-Z committed Oct 19, 2022
1 parent 704cf13 commit 16100dc
Showing 1 changed file with 150 additions and 61 deletions.
211 changes: 150 additions & 61 deletions MILP-demo.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,15 +9,25 @@
},
{
"cell_type": "code",
"execution_count": 22,
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"from scipy.optimize import milp, LinearConstraint\n",
"from qiskit_optimization import QuadraticProgram\n",
"from qiskit import QiskitError\n",
"import numpy as np\n",
"#from qiskit_optimization.algorithms import MinimumEigenOptimizer\n"
"from qiskit_optimization.algorithms import GurobiOptimizer\n"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
Expand All @@ -29,79 +39,133 @@
"### Example 1\n",
"\n",
"\n",
"$max \\ x_1$\n",
"$${\\rm max} \\ x_1+x_2+2x_3 \\\\\n",
"\n",
"$-x_1+x_2 \\leq 1 \\\\ 3x_1+2x_2 \\leq 12 \\\\ 2x_1+3x_2 \\leq 12$\n",
"-x_1+x_2+x_3 \\leq 1 \\\\ 3x_1+2x_2 - x_3 \\leq 12 \\\\ x_1+x_2 +x_3 \\leq 10 \\\\ 2x_1+3x_2 + 3x_3 \\geq 12 \\\\ \n",
"\n",
"$x_1, x_2 \\in \\mathcal{Z}$"
"x_1, x_2 \\in \\mathcal{Z}\n",
"$$\n",
"\n",
"Let's solve it with the milp solver!"
]
},
{
"cell_type": "code",
"execution_count": 4,
"execution_count": 14,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"x1: 1.0 x2: 2.0\n"
" fun: -15.0\n",
" message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'\n",
" mip_dual_bound: -15.0\n",
" mip_gap: 0.0\n",
" mip_node_count: 1\n",
" status: 0\n",
" success: True\n",
" x: array([5., 0., 5.])\n"
]
}
],
"source": [
"c = -np.array([0, 1])\n",
"A = np.array([[-1, 1], [3, 2], [2, 3]])\n",
"b_u = np.array([1, 12, 12])\n",
"b_l = np.full_like(b_u, -np.inf)\n",
"c = -np.array([1, 2, 2])\n",
"A = np.array([[-1, 1, 1], [3, 2, -1], [1, 1, 1], [2, 3, 3]])\n",
"b_u = np.array([1, 12, 10, np.inf])\n",
"b_l = np.array([-np.inf, -np.inf, -np.inf, 12])\n",
"\n",
"constraints = LinearConstraint(A, b_l, b_u)\n",
"integrality = np.ones_like(c)\n",
"integrality = np.array([1, 1, 0])\n",
"#integrality = np.zeros(3)\n",
"res = milp(c=c, constraints=constraints, integrality=integrality)\n",
"\n",
"print('x1:', res.x[0], 'x2: ', res.x[1])\n"
"print(res)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We see the MILP solver tells us the object function is maximized when $x_1 = 6, x_2 = 0, x_3 = 7$.\n",
"\n",
"As the next step, let's create a corresponding `QuadraticProgram` class in `qiskit_optimization`."
]
},
{
"cell_type": "code",
"execution_count": 41,
"execution_count": 17,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\\ This file has been generated by DOcplex\n",
"\\ ENCODING=ISO-8859-1\n",
"\\Problem name: example-1\n",
"Problem name: example-1\n",
"\n",
"Maximize\n",
" obj: x2\n",
"Subject To\n",
" c0: - x1 + x2 <= 1\n",
" c1: 3 x1 + 2 x2 <= 12\n",
" c2: 2 x1 + 3 x2 <= 12\n",
" x1 + x2 + 2*x3\n",
"\n",
"Subject to\n",
" Linear constraints (4)\n",
" -x1 + x2 + x3 <= 1 'c0'\n",
" x1 + x2 + x3 <= 10 'c1'\n",
" 3*x1 + 2*x2 - x3 <= 12 'c2'\n",
" 2*x1 + 3*x2 + 3*x3 >= 12 'c3'\n",
"\n",
"Bounds\n",
" Integer variables (2)\n",
" 0 <= x1\n",
" 0 <= x2\n",
"\n",
"Generals\n",
" x1 x2\n",
"End\n",
" Continuous variables (1)\n",
" 0 <= x3\n",
"\n"
]
}
],
"source": [
"## Consider the possibility of using sparse matrix(to_dict)\n",
"\n",
"qp = QuadraticProgram('example-1')\n",
"qp.integer_var(name='x1')\n",
"qp.integer_var(name='x2')\n",
"qp.maximize(linear={'x2': 1})\n",
"qp.continuous_var(name='x3')\n",
"#qp.binary_var_list(2000)\n",
"qp.maximize(linear={'x1': 1, 'x2': 1, 'x3': 2})\n",
"#qp.minimize(linear={'x2': 1})\n",
"qp.linear_constraint({'x1': -1, 'x2': 1},'<=',1)\n",
"qp.linear_constraint({'x1': 3, 'x2': 2},'<=',12)\n",
"qp.linear_constraint({'x1': 2, 'x2': 3},'<=',12)\n",
"qp.linear_constraint({'x1': -1, 'x2': 1, 'x3': 1},'<=',1)\n",
"qp.linear_constraint({'x1': 1, 'x2': 1, 'x3': 1},'<=',10)\n",
"qp.linear_constraint({'x1': 3, 'x2': 2, 'x3': -1},'<=',12)\n",
"qp.linear_constraint({'x1': 2, 'x2': 3, 'x3': 3},'>=',12)\n",
"\n",
"print(qp)"
"print(qp.prettyprint())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We solve the `qp(QuadraticProgram)` with the Gurobi solver."
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"objective function value: 15.0\n",
"variable values: x1=5.0, x2=-0.0, x3=5.0\n",
"status: SUCCESS\n"
]
}
],
"source": [
"gurobi_result = GurobiOptimizer().solve(qp)\n",
"print(gurobi_result.prettyprint())"
]
},
{
Expand All @@ -115,9 +179,24 @@
},
{
"cell_type": "code",
"execution_count": 31,
"execution_count": 34,
"metadata": {},
"outputs": [],
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" fun: -15.0\n",
" message: 'Optimization terminated successfully. (HiGHS Status 7: Optimal)'\n",
" mip_dual_bound: -15.0\n",
" mip_gap: 0.0\n",
" mip_node_count: 1\n",
" status: 0\n",
" success: True\n",
" x: array([5., 0., 5.])\n"
]
}
],
"source": [
"## First, Check if qp is linear\n",
"\n",
Expand All @@ -134,48 +213,56 @@
"sense = qp.objective.sense.value\n",
"\n",
"## cost function for milp solver\n",
"c = qp.objective.linear.to_array() * sense\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([0., 1.])"
]
},
"execution_count": 34,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"qp.objective.linear.to_array()"
"c_qp = qp.objective.linear.to_array() * sense\n",
"\n",
"## constraints for milp solver\n",
"for i, constraint in enumerate(qp.linear_constraints):\n",
"\n",
" constraint_array = constraint.linear.to_array()\n",
" constraint_sense = constraint.sense.value # 0 for leq, 1 for geq\n",
" constraint_value = constraint.rhs\n",
" if i==0:\n",
" A_qp = constraint_array\n",
" bl_qp = np.array([constraint_value if constraint_sense==1 else -np.inf])\n",
" bu_qp = np.array([constraint_value if constraint_sense==0 else np.inf])\n",
" else:\n",
" A_qp = np.vstack((A_qp,constraint_array))\n",
" bl_qp = np.append(bl_qp, [constraint_value if constraint_sense==1 else -np.inf])\n",
" bu_qp = np.append(bu_qp, [constraint_value if constraint_sense==0 else np.inf])\n",
"\n",
"constraints = LinearConstraint(A_qp, bl_qp, bu_qp)\n",
"\n",
"## integrity for milp solver\n",
"integrality = np.array([ 1 if variable.vartype.value==2 else 0 for variable in qp.variables]) ## check on other vartype!\n",
"\n",
"raw_res = milp(c=c_qp, constraints=constraints, integrality=integrality)\n",
"\n",
"print(raw_res)\n",
" "
]
},
{
"cell_type": "code",
"execution_count": 60,
"execution_count": 29,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"array([-1., 1.])"
"array([[-1., 1., 1.],\n",
" [ 3., 2., -1.]])"
]
},
"execution_count": 60,
"execution_count": 29,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"qp.get_linear_constraint(0).linear.to_array()"
"A = qp.get_linear_constraint(0).linear.to_array()\n",
"B = qp.get_linear_constraint(1).linear.to_array()\n",
"A = np.vstack((A, B))\n",
"A"
]
},
{
Expand All @@ -200,10 +287,12 @@
},
{
"cell_type": "code",
"execution_count": null,
"execution_count": 26,
"metadata": {},
"outputs": [],
"source": []
"source": [
"integrality = np.array([ 1 if variable.vartype.value==2 else 0 for variable in qp.variables])"
]
}
],
"metadata": {
Expand Down

0 comments on commit 16100dc

Please sign in to comment.