1
$\begingroup$

enter image description here Here is my optimization problem I solve with classical optimization using DocPlex but I want to solve using quadratic program using Qiskit and QAOA as solver optimumeigen.

The requirement is that I need to solve it with quantum algorithms. For solving QUBO problems, I could use QAOA. But for solving quadratic programming problems with only continuous variables, I don't know which quantum algorithm I could use. Or maybe this problem could be converted to a QUBO problem and then handled by QAOA. I'm not sure.

enter image description here

$\endgroup$
1
  • 1
    $\begingroup$ Welcome to QCSE! Please add your code as text instead of using an image! $\endgroup$
    – FDGod
    Commented Jan 25 at 17:28

1 Answer 1

3
$\begingroup$

QUBO to QAOA/ADMM

To solve a continuous variable problem like yours, you should use the ADMM optimizer, and here is how you should do it.

1. Keep your code in making the classical cplex model as it is, i.e

import numpy as np
from docplex.mp.model import Model

ts=24 # simulation time

deltat = 1
Rgh = 5
Cgh = 6
uj = 3 #Given parameter

Tou_h = np.array([37.0, 38.5, 37.6, 36.6, 38.3, 37.9, 37.5, 38.6, 38.1, 37.3, 35.0, 36.3, 37.9, 36.4, 38.3, 37.7, 36.8, 35.8, 35.9, 36.0, 37.4, 37.3, 37.2, 37.5])
Tin_ref = np.array([26] * 10 + [24] * 14)

# Create a Docplex model
model = Model (name='temperature_minimization')

# Define decision variables
Tin_h = {t: model.continuous_var(name=f'Tin_h_{t}', lb=0) for t in range(ts)}
Uh_gh = {t: model.continuous_var (name=f'Uh_gh_{t}', lb=-15.6, ub=15.6) for t in range(ts)}

# Define constraints
for t in range(ts - 1):
  model.add_constraint (Tin_h[t + 1] ==Tin_h[t] - (deltat / Rgh * Cgh) * (Tin_h[t] - Tou_h[t]) - uj * Rgh * Uh_gh[t])



# Define an objective function
obj_function = model.sum((Tin_h[t] - Tin_ref[t]) ** 2 for t in range(ts))
model.minimize(obj_function)

# solving the model
model.solve()

# Display the results.
print("Objective Value:", model.objective_value)
for t in range(ts):
  print (f'Tin_h_{t} = {Tin_h[t].solution_value}, Uh_gh_{t} = {Uh_gh[t].solution_value}')

2. Make sure to install required qiskit dependencies using:

%pip install qiskit
%pip install cplex
%pip install docplex
%pip install qiskit_optimization

and then make the necessary imports for the qiskit package:

from qiskit_optimization import QuadraticProgram
from qiskit_optimization.translators.docplex_mp import to_docplex_mp
import matplotlib.pyplot as plt

from docplex.mp.model import Model

from qiskit_algorithms import QAOA, NumPyMinimumEigensolver
from qiskit_algorithms.optimizers import COBYLA
from qiskit.primitives import Sampler
from qiskit_optimization.algorithms import CobylaOptimizer, MinimumEigenOptimizer
from qiskit_optimization.algorithms.admm_optimizer import ADMMParameters, ADMMOptimizer
from qiskit_optimization.translators import from_docplex_mp

Do note that this comes with a free version of cplex limited to solving $1000$ variables and $1000$ constraints. If your problem exceeds that, it won't solve it.

3. Import your model such that qiskit can understand it:

# load from a Docplex model
mod = from_docplex_mp(model)
print(type(mod))
print()
print(mod.prettyprint())

you can check the details about your model like this, it displays all the linear variables, constants, the cost matrix, everything:

print("constant:\t\t\t", mod.objective.constant)
print("linear dict:\t\t\t", mod.objective.linear.to_dict())
print("linear array:\t\t\t", mod.objective.linear.to_array())
print("linear array as sparse matrix:\n", mod.objective.linear.coefficients, "\n")
print("quadratic dict w/ index:\t", mod.objective.quadratic.to_dict())
print("quadratic dict w/ name:\t\t", mod.objective.quadratic.to_dict(use_name=True))
print(
    "symmetric quadratic dict w/ name:\t",
    mod.objective.quadratic.to_dict(use_name=True, symmetric=True),
)
print("quadratic matrix:\n", mod.objective.quadratic.to_array(), "\n")
print("symmetric quadratic matrix:\n", mod.objective.quadratic.to_array(symmetric=True), "\n")
print("quadratic matrix as sparse matrix:\n", mod.objective.quadratic.coefficients)

4. Define the solvers you'll make use in solving the problem:

Classical Solution

3-ADMM-H needs a QUBO optimizer to solve the QUBO subproblem, and a continuous optimizer to solve the continuous convex constrained subproblem. We first solve the problem classically:

we use the MinimumEigenOptimizer with the NumPyMinimumEigenSolver as a classical and exact QUBO solver and we use the CobylaOptimizer as a continuous convex solver. 3-ADMM-H supports any other suitable solver available in Qiskit optimization.

# define COBYLA optimizer to handle convex continuous problems.
cobyla = CobylaOptimizer()

# define QAOA via the minimum eigen optimizer
qaoa = MinimumEigenOptimizer(QAOA(sampler=Sampler(), optimizer=COBYLA()))

# exact QUBO solver as a classical benchmark
exact = MinimumEigenOptimizer(NumPyMinimumEigensolver())  # to solve QUBOs

# In case CPLEX is installed it can also be used for the convex problems, the QUBO,
# or as a benchmark for the full problem.
#
# cplex = CplexOptimizer()

Parameters

The 3-ADMM-H are wrapped in class ADMMParameters. Customized parameter values can be set as arguments of the class.

admm_params = ADMMParameters(
    rho_initial=1001, beta=1000, factor_c=900, maxiter=100, three_block=True, tol=1.0e-6
)

Calling 3-ADMM-H algorithm

To invoke the 3-ADMM-H algorithm, an instance of the ADMMOptimizer class needs to be created. This takes ADMM-specific parameters and the subproblem optimizers separately into the constructor. The solution returned is an instance of the OptimizationResult class.

# define QUBO optimizer
qubo_optimizer = exact
# qubo_optimizer = cplex  # uncomment to use CPLEX instead

# define classical optimizer
convex_optimizer = cobyla
# convex_optimizer = cplex  # uncomment to use CPLEX instead

# initialize ADMM with classical QUBO and convex optimizer
admm = ADMMOptimizer(
    params=admm_params, qubo_optimizer=qubo_optimizer, continuous_optimizer=convex_optimizer
)
  1. Run the ADMM Solver for your problem
# run ADMM to solve problem
result = admm.solve(mod)

Classical Solver Result

The 3-ADMM-H solution can be then printed and visualized. The x attribute of the solution contains respectively, the values of the binary decision variables and the values of the continuous decision variables. The fval is the objective value of the solution.

print(result.prettyprint())

It will give an output like this:

objective function value: 1354.8542360685278
variable values: Tin_h_0=17.663029761048875, Tin_h_1=18.491408258771582, Tin_h_2=17.54786970153851, Tin_h_3=18.6364987479822, Tin_h_4=17.496378700095512, Tin_h_5=17.889667346808253, Tin_h_6=18.780309372587702, Tin_h_7=18.77321307560003, Tin_h_8=19.053101875660403, Tin_h_9=17.679025703596793, Tin_h_10=17.291236701341333, Tin_h_11=15.782773277770197, Tin_h_12=16.724625174535355, Tin_h_13=16.049729887057467, Tin_h_14=17.336396869454838, Tin_h_15=16.592380862605346, Tin_h_16=17.41354648393138, Tin_h_17=15.56175113614779, Tin_h_18=16.207082861820027, Tin_h_19=16.577530722487364, Tin_h_20=16.749629149997062, Tin_h_21=17.634392099010093, Tin_h_22=16.643925784908333, Tin_h_23=17.850220826905257, Uh_gh_0=1.4917323859345744, Uh_gh_1=1.6635899097804796, Uh_gh_2=1.5315951541140045, Uh_gh_3=1.5130881033538714, Uh_gh_4=1.6380704608781755, Uh_gh_5=1.5414504772033761, Uh_gh_6=1.4980483366588286, Uh_gh_7=1.5674837006146383, Uh_gh_8=1.6153569280847424, Uh_gh_9=1.5955305438626206, Uh_gh_10=1.51726529213077, Uh_gh_11=1.578588011327375, Uh_gh_12=1.739023005202363, Uh_gh_13=1.5422438102089115, Uh_gh_14=1.7266893175669118, Uh_gh_15=1.633865156236503, Uh_gh_16=1.674369304471061, Uh_gh_17=1.5760377940633585, Uh_gh_18=1.5507368470099077, Uh_gh_19=1.5423243137003648, Uh_gh_20=1.5930454713993645, Uh_gh_21=1.6392797196859752, Uh_gh_22=1.564066267740873, Uh_gh_23=2.1134932838436886
status: SUCCESS

6. Quantum Solution

We now solve the same optimization problem with QAOA as a QUBO optimizer, running on the simulated quantum device. First, one needs to select the classical optimizer of the eigensolver QAOA. Then, the simulation backend is set. Finally, the eigensolver is wrapped into the MinimumEigenOptimizer class. A new instance of ADMMOptimizer is populated with QAOA as QUBO optimizer.

# define QUBO optimizer
qubo_optimizer = qaoa

# define classical optimizer
convex_optimizer = cobyla
# convex_optimizer = cplex  # uncomment to use CPLEX instead

# initialize ADMM with quantum QUBO optimizer and classical convex optimizer
admm_q = ADMMOptimizer(
    params=admm_params, qubo_optimizer=qubo_optimizer, continuous_optimizer=convex_optimizer
)

then run ADMM to solve it

# run ADMM to solve problem
result_q = admm_q.solve(mod)

Quantum Solver Results

Here we present the results obtained from the quantum solver. As in the example above x stands for the solution, and the fval is for objective value.

print(result.prettyprint())

You'll get an output like this:

objective function value: 1354.8542360685278
variable values: Tin_h_0=17.663029761048875, Tin_h_1=18.491408258771582, Tin_h_2=17.54786970153851, Tin_h_3=18.6364987479822, Tin_h_4=17.496378700095512, Tin_h_5=17.889667346808253, Tin_h_6=18.780309372587702, Tin_h_7=18.77321307560003, Tin_h_8=19.053101875660403, Tin_h_9=17.679025703596793, Tin_h_10=17.291236701341333, Tin_h_11=15.782773277770197, Tin_h_12=16.724625174535355, Tin_h_13=16.049729887057467, Tin_h_14=17.336396869454838, Tin_h_15=16.592380862605346, Tin_h_16=17.41354648393138, Tin_h_17=15.56175113614779, Tin_h_18=16.207082861820027, Tin_h_19=16.577530722487364, Tin_h_20=16.749629149997062, Tin_h_21=17.634392099010093, Tin_h_22=16.643925784908333, Tin_h_23=17.850220826905257, Uh_gh_0=1.4917323859345744, Uh_gh_1=1.6635899097804796, Uh_gh_2=1.5315951541140045, Uh_gh_3=1.5130881033538714, Uh_gh_4=1.6380704608781755, Uh_gh_5=1.5414504772033761, Uh_gh_6=1.4980483366588286, Uh_gh_7=1.5674837006146383, Uh_gh_8=1.6153569280847424, Uh_gh_9=1.5955305438626206, Uh_gh_10=1.51726529213077, Uh_gh_11=1.578588011327375, Uh_gh_12=1.739023005202363, Uh_gh_13=1.5422438102089115, Uh_gh_14=1.7266893175669118, Uh_gh_15=1.633865156236503, Uh_gh_16=1.674369304471061, Uh_gh_17=1.5760377940633585, Uh_gh_18=1.5507368470099077, Uh_gh_19=1.5423243137003648, Uh_gh_20=1.5930454713993645, Uh_gh_21=1.6392797196859752, Uh_gh_22=1.564066267740873, Uh_gh_23=2.1134932838436886
status: SUCCESS
$\endgroup$
2
  • $\begingroup$ Thank you so much for your help, but I don't know why the optimization using your code doesn't give a good result like a classical optimization using Pyomo and Ipopt as solvers. and if possible, I need to solve this problem using QAOA and VQE. not Cobyla as a classical solver that you use in part of the quantum solution. $\endgroup$
    – said zhn
    Commented Jan 29 at 11:45
  • $\begingroup$ This is QAOA, and QAOA uses classical optimizer to solve the problem. The classical and quantum results are exactly the same as you can see from the code. $\endgroup$ Commented Feb 1 at 15:07

Not the answer you're looking for? Browse other questions tagged or ask your own question.