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(line_seperator)
print("N(Number of assets): ", N)
print("M(Number of levels): ", M)
print(line_seperator)
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(line_seperator)
print("R(Risk tolerance parameter): ", R)
print("B(Total budget): ", B)
print("D(Diversification requirement): ", D)
print("K(Cardinality constraint): ", K)
print(line_seperator)
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)
print(qp.export_as_lp_string())
# General imports
import numpy as np
import warnings
warnings.filterwarnings("ignore")
# 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)
backend.name
op,offset = qp.to_ising()
print("offset: {}".format(offset))
print("operator:")
print(op)
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)
hamiltonian_isa
def cost_func(params, ansatz, hamiltonian, estimator):
"""Return estimate of energy from estimator
Parameters:
params (ndarray): Array of ansatz parameters
ansatz (QuantumCircuit): Parameterized ansatz circuit
hamiltonian (SparsePauliOp): Operator representation of Hamiltonian
estimator (EstimatorV2): Estimator primitive instance
Returns:
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.
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.measure_all()
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
session.close()
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(
channel='ibm_quantum',
instance='ibm-q/open/main',
token=''
)
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.
