As suggested in TranspilerError: 'Number of qubits (40) in QAOA is greater than maximum (30) in the coupling_map' issue. I managed to submit the job to IBM quantum and the circuit runs successfully.

The code is as follows.

import numpy as np
from docplex.mp.advmodel import AdvModel
from qiskit_optimization.translators import from_docplex_mp

mdl = AdvModel(name="Portfolio optimization")
N = 10   # Number of assets
M = 4   # Number of levels

lambda_risk = 5
lambda_linking = 10
lambda_budget = 20
lambda_diversification = 15
lambda_cardinality = 10

R = 100   # Risk tolerance parameter
B = 50    # Total budget
D = 16    # Diversification requirement
K = 6     # Cardinality constraint

# Vector of mu_i values, expected returns for each asset 
mu = np.random.uniform(low=0.1, high=1.0, size=N)

# sigma_ij values, representing the covariance between assets i and j
sigma = np.random.uniform(low=0.01, high=0.1, size=(N, N))
# c_i values, representing the cost of each asset
c = np.random.uniform(low=1, high=10, size=N)
line_seperator = "----------------------------------------"
print("N(Number of assets): ", N)
print("M(Number of levels): ", M)
print("lambda_risk(Risk penalty coefficient): ", lambda_risk)
print("lambda_linking(Linking penalty coefficient): ", lambda_linking)
print("lambda_budget(Budget penalty coefficient): ", lambda_budget)
print("lambda_diversification(Diversification penalty coefficient): ", lambda_diversification)
print("lambda_cardinality(Cardinality penalty coefficient): ", lambda_cardinality)
print("R(Risk tolerance parameter): ", R)
print("B(Total budget): ", B)
print("D(Diversification requirement): ", D)
print("K(Cardinality constraint): ", K)
print("mu(Expected returns for each asset): ", mu)
print("sigma(Covariance between assets): ", sigma)
print("c(Cost of each asset): ", c)

z = {(i, l): mdl.binary_var(name=f"z_{i}_{l}") for i in range(N) for l in range(M)}
y = [mdl.binary_var(name=f"y_{i}") for i in range(N)]

#objective function
obj = mdl.sum(mu[i] * mdl.sum(2**l * z[i, l] for l in range(M)) for i in range(N))

penalty_risk = mdl.sum(mdl.sum(2**l * z[i, l] for l in range(M)) * sigma[i][j] * mdl.sum(2**l * z[j, l] for l in range(M)) for i in range(N) for j in range(N)) - R
# risk_term = lambda_risk * penalty_risk ** 2
risk_term = lambda_risk * penalty_risk

penalty_budget = mdl.sum(c[i] * mdl.sum(2**l * z[i, l] for l in range(M)) for i in range(N)) - B
# budget_term = lambda_budget * penalty_budget ** 2
budget_term = lambda_budget * penalty_budget

penalty_diversification = mdl.sum(mdl.sum(2**l * z[i, l] for l in range(M)) for i in range(N)) - D
# diversification_term = lambda_diversification * penalty_diversification ** 2
diversification_term = lambda_diversification * penalty_diversification

penalty_cardinality = mdl.sum(y[i] for i in range(N)) - K
# cardinality_term = lambda_cardinality * penalty_cardinality ** 2
cardinality_term = lambda_cardinality * penalty_cardinality

penalty_linking = mdl.sum(mdl.sum(2**l * z[i, l] for l in range(M)) for i in range(N)) - D * mdl.sum(y[i] for i in range(N))
# linking_term = lambda_linking * penalty_linking ** 2
linking_term = lambda_linking * penalty_linking

mdl.minimize(risk_term + budget_term+ diversification_term + cardinality_term + linking_term - obj)
qp = from_docplex_mp(mdl)
# General imports
import numpy as np
import warnings


# Pre-defined ansatz circuit, operator class and visualization tools
from qiskit.quantum_info import SparsePauliOp
from qiskit.visualization import plot_distribution

# Qiskit Runtime
from qiskit_ibm_runtime import QiskitRuntimeService, Session
from qiskit_ibm_runtime import EstimatorV2 as Estimator
from qiskit_ibm_runtime import SamplerV2 as Sampler

# SciPy minimizer routine
from scipy.optimize import minimize

# rustworkx graph library
import rustworkx as rx
from rustworkx.visualization import mpl_draw

# To run on hardware, select the backend with the fewest number of jobs in the queue
service = QiskitRuntimeService(channel="ibm_quantum",token="Token")
backend = service.least_busy(min_num_qubits=40,operational=True, simulator=False)
op,offset = qp.to_ising()
print("offset: {}".format(offset))

from qiskit.circuit.library import QAOAAnsatz
ansatz = QAOAAnsatz(op, reps=2)

ansatz.decompose(reps=3).draw(output="mpl", style="iqp")

from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager

target = backend.target
pm = generate_preset_pass_manager(target=target, optimization_level=3)

ansatz_isa = pm.run(ansatz)

hamiltonian_isa = op.apply_layout(ansatz_isa.layout)
def cost_func(params, ansatz, hamiltonian, estimator):
    """Return estimate of energy from estimator

        params (ndarray): Array of ansatz parameters
        ansatz (QuantumCircuit): Parameterized ansatz circuit
        hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
        estimator (EstimatorV2): Estimator primitive instance

        float: Energy estimate
    pub = (ansatz, [hamiltonian], [params])
    result = estimator.run(pubs=[pub]).result()
    cost = result[0].data.evs[0]

    return cost
session = Session(backend=backend)

# Configure estimator
estimator = Estimator(session=session)
estimator.options.default_shots = 10_000
estimator.options.dynamical_decoupling.enable = True

# Configure sampler
sampler = Sampler(session=session)
sampler.options.default_shots = 10_000
sampler.options.dynamical_decoupling.enable = True

x0 = 2 * np.pi * np.random.rand(ansatz_isa.num_parameters)
res = minimize(cost_func, x0, args=(ansatz_isa, hamiltonian_isa, estimator), method="COBYLA")

But the Python interpreter runs into an error due to the long time it takes to solve this. But when I go to IBM Quantum platform the job has finished.enter image description here

The suggested way to move forward acoording to IBM QAOA is in the link as follows.

# Assign solution parameters to ansatz
qc = ansatz.assign_parameters(res.x)
# Add measurements to our circuit
qc_isa = pm.run(qc)
qc_isa.draw(output="mpl", idle_wires=False, style="iqp")

result = sampler.run([qc_isa]).result()
samp_dist = result[0].data.meas.get_counts()
# Close the session since we are now done with it

But since my interpretter times out I don't get to use the res.x parameters.

But the IBM platforms provides a way to use the results of the job as follows.

from qiskit_ibm_runtime import QiskitRuntimeService

service = QiskitRuntimeService(
job = service.job('cqzvw6gktf3g00884pz0')
job_result = job.result()

for idx, pub_result in enumerate(job_result):
    print(f"Expectation values for pub {idx}: {pub_result.data.evs}")

I want to know how I can use this results to carry forward with the calculation. I am new to QAOA so this might be straight forward. Any help is valued.



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.