1
$\begingroup$

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.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.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.

$\endgroup$

0

Your Answer

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