# Qiskit Optimizer 範例: 二次函數求最佳解
Adapted from [IBM Quantum Challenge Africa 2021](https://challenges.quantum-computing.ibm.com/africa21) Lab 1
中文敘述部分：2021 [(CC BY-NC-SA 4.0)](https://creativecommons.org/licenses/by-nc-sa/4.0/deed.en) Elton Huang
敘述為強調流程的重點，因而簡化了許多，建議大家還是要一起看原文檔
這個練習示範高階的量子計算使用模式：如何使用 `Qiskit Optimization` 應用模組對二次函數問題求最佳解，包括
1. 如何以 Qiskit 定義傳統與數學問題,
2. 如何定義量子算法,
3. 如何針對所定義的問題執行量子算法解題,
4. 如何以模擬和雲端實機的方式執行算法解題.
```python
# 先匯入需用的模組
# Import auxiliary libraries
import numpy as np
# Import Qiskit
from qiskit import IBMQ, Aer
from qiskit.algorithms import QAOA, VQE, NumPyMinimumEigensolver
from qiskit.algorithms.optimizers import COBYLA
from qiskit.utils import QuantumInstance, algorithm_globals
from qiskit.providers.aer.noise.noise_model import NoiseModel
from qiskit_optimization import QuadraticProgram
from qiskit_optimization.algorithms import MinimumEigenOptimizer
from qiskit_optimization.converters import QuadraticProgramToQubo
import qiskit.test.mock as Fake
```
## 問題敘述與建模
有一款單機遊戲，對手是電腦，得分完全取決與你組合中角色 (包括主教 W、巫師 S、騎士 M、士兵 P) 的配置：
1. 主教 W 得 2 分
2. 巫師 S 得 1 分
3. 騎士 M 得 4 分
4. 一位主教 W 搭配一位巫師 S 得 2.4 分
5. 一位主教 W 搭配一位騎士 M 得 4 分
6. 一位主教 W 搭配一位士兵 P 得 4 分
7. 一位巫師 S 搭配一位騎士 M 得 2 分
8. 一位巫師 S 搭配一位士兵 P 得 1 分
9. 一位騎士 M 搭配一位士兵 P 得 5 分
你可以安排的組合有以下限制
1. 總共不能超過 3 個人
2. 每種角色不能超過 1 位
對於這類的問題，我們可以用一個二次函數的模型來描述，之後就可以操作程式的工具的解題。如果得分是 $f$:
$f = 2w + s + 4m + 2.4ws + 4wm + 4wp + 2sm + sp + 5mp$
$w + s + m + p \le 3$
$0 \le w \le 1$
$0 \le s \le 1$
$0 \le m \le 1$
$0 \le p \le 1$
如果用 Qikit Optimization 套件的工具來描述，則如下:
(程式中變數名稱請取第一個字母來對應題敘)
```python
def cropyield_quadratic_program():
cropyield = QuadraticProgram(name="Crop Yield") # 先給這個問題個名稱
##############################
# Put your implementation here
cropyield.integer_var(name="Wheat", lowerbound=0, upperbound=1)
cropyield.integer_var(name="Soybeans", lowerbound=0, upperbound=1)
cropyield.integer_var(name="Maize", lowerbound=0, upperbound=1)
cropyield.integer_var(name="PushPull", lowerbound=0, upperbound=1)
cropyield.maximize(
linear={"Wheat": 2, "Soybeans": 1, "Maize": 4},
quadratic={("Wheat", "Soybeans"): 2.4, ("Wheat", "Maize"): 4, ("Wheat", "PushPull"): 4, ("Soybeans", "Maize"): 2, ("Soybeans", "PushPull"): 1, ("Maize", "PushPull"): 5},
)
cropyield.linear_constraint(linear={"Wheat": 1, "Soybeans": 1, "Maize": 1, "PushPull": 1}, sense="<=", rhs=3)
#
##############################
return cropyield
cropyield = cropyield_quadratic_program()
```
下面的操作將我們用 Qiskit 建立的模型整理出來，可以用來確認模型沒有錯誤
```python
print(cropyield.export_as_lp_string())
```
\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: Crop Yield
Maximize
obj: 2 Wheat + Soybeans + 4 Maize + [ 4.800000000000 Wheat*Soybeans
+ 8 Wheat*Maize + 8 Wheat*PushPull + 4 Soybeans*Maize
+ 2 Soybeans*PushPull + 10 Maize*PushPull ]/2
Subject To
c0: Wheat + Soybeans + Maize + PushPull <= 3
Bounds
Wheat <= 1
Soybeans <= 1
Maize <= 1
PushPull <= 1
Generals
Wheat Soybeans Maize PushPull
End
下面的操作估計我們需要幾個 qubits 來求解
```python
# Estimate the number of qubits required
ising_operations, _ = (
QuadraticProgramToQubo()
.convert(
cropyield,
)
.to_ising()
)
print(f"Number of qubits required is {ising_operations.num_qubits}")
```
Number of qubits required is 6
最後再做一個轉換將變數轉成二元變數 (細節可以之後再了解)
```python
QuadraticProgramToQubo().convert(cropyield)
```
\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: Crop Yield
Minimize
obj: - 160.400000000000 Wheat@0 - 159.400000000000 Soybeans@0
- 162.400000000000 Maize@0 - 158.400000000000 PushPull@0
- 158.400000000000 c0@int_slack@0 - 316.800000000000 c0@int_slack@1 + [
52.800000000000 Wheat@0^2 + 100.800000000000 Wheat@0*Soybeans@0
+ 97.600000000000 Wheat@0*Maize@0 + 97.600000000000 Wheat@0*PushPull@0
+ 105.600000000000 Wheat@0*c0@int_slack@0
+ 211.200000000000 Wheat@0*c0@int_slack@1 + 52.800000000000 Soybeans@0^2
+ 101.600000000000 Soybeans@0*Maize@0
+ 103.600000000000 Soybeans@0*PushPull@0
+ 105.600000000000 Soybeans@0*c0@int_slack@0
+ 211.200000000000 Soybeans@0*c0@int_slack@1 + 52.800000000000 Maize@0^2
+ 95.600000000000 Maize@0*PushPull@0
+ 105.600000000000 Maize@0*c0@int_slack@0
+ 211.200000000000 Maize@0*c0@int_slack@1 + 52.800000000000 PushPull@0^2
+ 105.600000000000 PushPull@0*c0@int_slack@0
+ 211.200000000000 PushPull@0*c0@int_slack@1
+ 52.800000000000 c0@int_slack@0^2
+ 211.200000000000 c0@int_slack@0*c0@int_slack@1
+ 211.200000000000 c0@int_slack@1^2 ]/2 + 237.600000000000
Subject To
Bounds
0 <= Wheat@0 <= 1
0 <= Soybeans@0 <= 1
0 <= Maize@0 <= 1
0 <= PushPull@0 <= 1
0 <= c0@int_slack@0 <= 1
0 <= c0@int_slack@1 <= 1
Binaries
Wheat@0 Soybeans@0 Maize@0 PushPull@0 c0@int_slack@0 c0@int_slack@1
End
## 本機模擬解題
模型建好後，有以下三種方式用 Qikit 程式工具解題：
1. 使用本機模擬器
2. 使用 IBM 雲端模擬器
3. 使用 IBM Quantum 量子電腦
這三種都稱為 backend (後端)，都可以用來解題
### 傳統的解法
不過使用量子計算解題之前，我們先看看用傳統的計算模式解題的結果
```python
def get_classical_solution_for(quadprog: QuadraticProgram):
# Create solver
solver = NumPyMinimumEigensolver()
# Create optimizer for solver
optimizer = MinimumEigenOptimizer(solver)
# Return result from optimizer
return optimizer.solve(quadprog)
# Get classical result
classical_result = get_classical_solution_for(cropyield)
# Format and print result
print("Solution found using the classical method:\n")
print(f"Maximum crop-yield is {classical_result.fval} tons")
print(f"Crops used are: ")
_crops = [v.name for v in cropyield.variables]
for cropIndex, cropHectares in enumerate(classical_result.x):
print(f"\t{cropHectares} ha of {_crops[cropIndex]}")
```
Solution found using the classical method:
Maximum crop-yield is 19.0 tons
Crops used are:
1.0 ha of Wheat
0.0 ha of Soybeans
1.0 ha of Maize
1.0 ha of PushPull
### 使用本機模擬器解題
接著用量子計算來解題。第一種方式：使用本機模擬器。
```python
# We will use the Aer provided QASM simulator
backend = Aer.get_backend("qasm_simulator")
# Given we are using a simulator, we will fix the algorithm seed to ensure our results are reproducible
algorithm_globals.random_seed = 271828
```
#### 量子計算 QAOA 解題 (本機模擬)
接下來有分幾種模式，先用 QAOA
```python
def get_QAOA_solution_for(
quadprog: QuadraticProgram, quantumInstance: QuantumInstance, optimizer=None,
):
_eval_count = 0
def callback(eval_count, parameters, mean, std):
nonlocal _eval_count
_eval_count = eval_count
# Create solver
solver = QAOA(
optimizer=optimizer, quantum_instance=quantumInstance, callback=callback,
)
# Create optimizer for solver
optimizer = MinimumEigenOptimizer(solver)
# Get result from optimizer
result = optimizer.solve(quadprog)
return result, _eval_count
# Create a QuantumInstance
simulator_instance = QuantumInstance(
backend=backend,
seed_simulator=algorithm_globals.random_seed,
seed_transpiler=algorithm_globals.random_seed,
)
# Get QAOA result
qaoa_result, qaoa_eval_count = get_QAOA_solution_for(cropyield, simulator_instance)
# Format and print result
print("Solution found using the QAOA method:\n")
print(f"Maximum crop-yield is {qaoa_result.fval} tons")
print(f"Crops used are: ")
for cropHectares, cropName in zip(qaoa_result.x, qaoa_result.variable_names):
print(f"\t{cropHectares} ha of {cropName}")
print(f"\nThe solution was found within {qaoa_eval_count} evaluations of QAOA.")
```
Solution found using the QAOA method:
Maximum crop-yield is 19.0 tons
Crops used are:
1.0 ha of Wheat
0.0 ha of Soybeans
1.0 ha of Maize
1.0 ha of PushPull
The solution was found within 3 evaluations of QAOA.
結果和我們用傳統模式的方法所得到的是一樣的
#### 量子計算 VQE 解題 (本機模擬)
量子計算中第二種模式 VQE:
```python
def get_VQE_solution_for(
quadprog: QuadraticProgram, quantumInstance: QuantumInstance, optimizer=None,
):
_eval_count = 0
def callback(eval_count, parameters, mean, std):
nonlocal _eval_count
_eval_count = eval_count
# Create solver and optimizer
solver = VQE(
optimizer=optimizer, quantum_instance=quantumInstance, callback=callback
)
# Create optimizer for solver
optimizer = MinimumEigenOptimizer(solver)
# Get result from optimizer
result = optimizer.solve(quadprog)
return result, _eval_count
# Create a QuantumInstance
simulator_instance = QuantumInstance(
backend=backend,
seed_simulator=algorithm_globals.random_seed,
seed_transpiler=algorithm_globals.random_seed,
)
# Get VQE result
vqe_result, vqe_eval_count = get_VQE_solution_for(cropyield, simulator_instance)
# Format and print result
print("Solution found using the VQE method:\n")
print(f"Maximum crop-yield is {vqe_result.fval} tons")
print(f"Crops used are: ")
for cropHectares, cropName in zip(vqe_result.x, vqe_result.variable_names):
print(f"\t{cropHectares} ha of {cropName}")
print(f"\nThe solution was found within {vqe_eval_count} evaluations of VQE")
```
Solution found using the VQE method:
Maximum crop-yield is 19.0 tons
Crops used are:
1.0 ha of Wheat
0.0 ha of Soybeans
1.0 ha of Maize
1.0 ha of PushPull
The solution was found within 25 evaluations of VQE
結果和之前兩個方法所得到的是一樣的
(... 這裡略過模擬真實量子電腦 ...)
## 雲端執行
最上面設定問題時我們說每個角色最多 1 人，全部角色最多 3 人，現在我們把限制擴大到每個角色最多 10 人，全部角色最多也是 10 人，這樣需要幾個 qubits 呢？
```python
# Function to estimate the number of qubits required
def estimate_number_of_qubits_required_for(max_hectares_per_crop, hectares_available):
return int(
4 * np.ceil(np.log2(max_hectares_per_crop + 1))
+ np.ceil(np.log2(hectares_available + 1))
)
# Our new problem parameters
hectares_available = 10
max_hectares_per_crop = 10
# Retrieving the number of qubits required
number_of_qubits_required = estimate_number_of_qubits_required_for(
max_hectares_per_crop=max_hectares_per_crop, hectares_available=hectares_available
)
print(
f"Optimizing a {hectares_available} ha farm with each crop taking up to {max_hectares_per_crop} ha each,",
f"the computation is estimated to require {number_of_qubits_required} qubits.",
)
```
Optimizing a 10 ha farm with each crop taking up to 10 ha each, the computation is estimated to require 20 qubits.
The number of qubits required is related to the constraints in the quadratic program and how the integer variables are converted to binary variables. In fact, the scaling of the number of qubits, as a function of the hectares available, is logarithmic in nature; owing to this conversion.
接著我們看看如何透過雲端試著用實機來執行
(編按：我看得到雲端實機最多 7 quibts 只有一台 ibmq_perth，最多是 5 quibits，解題加上誤差，連最初的限制都不太夠用)
To use the IBM Quantum platform is easy. First you need to load the account you enabled in the week 0 content. If you didn't complete this, follow this [quick guide](https://quantum-computing.ibm.com/lab/docs/iql/manage/account/ibmq) on connecting your IBM Quantum account with Qiskit in python and Jupyter.
```python
IBMQ.load_account()
```
ibmqfactory.load_account:WARNING:2021-09-18 05:56:19,374: Credentials are already in use. The existing account in the session will be replaced.
<AccountProvider for IBMQ(hub='ibm-q', group='open', project='main')>
IBM Quantum backends are accessed through a provider, which manages the devices to which you have access. For this challenge, you have access to the new `ibm_perth` quantum computer! Typically, you would find your provider details under your [IBM Quantum account details](https://quantum-computing.ibm.com/account). Under your account you can see the different hubs, groups, and projects you are a part of. Qiskit allows us to retrieve a provider using just the hub, group, and project as follows:
```python
provider = IBMQ.get_provider(hub="ibm-q", group="open", project="main")
```
However, because we have given you special access for this challenge, we are going to retrive the provider using a different method. Execute the code cell below to retrieve the correct provider.
```python
provider = None
for prov in IBMQ.providers():
if (
"iqc-africa-21" in prov.credentials.hub
and "q-challenge" in prov.credentials.group
and "ex1" in prov.credentials.project
):
# Correct provider found
provider = prov
if provider == None:
print("ERROR: The expected provider was not found!")
else:
print("Yay! The expected provider was found!")
```
Yay! The expected provider was found!
If the above code cell returned an error, you may not yet have access to the real quantum computer. The list of participants is updated daily, so you may have to wait some time before the correct provider appears. If you need assistance, send a message to the challenge Slack channel [#challenge-africa-2021](https://qiskit.slack.com/archives/C02C8MKP153) and make sure to tag the admin team with [@africa_admin](#).
------
To retrieve a backend from the provider, one needs only request it by name. For example, we can request `ibm_perth` as follows.
```python
backend_real = provider.get_backend("ibm_perth")
```
(... 省略關於雲端實機與模擬選擇的敘述 ...)
We can create a new `QuantumInstance` object to contain our real quantum computer backend, similar to how we created one to manage our simulator. For real devices, there is an extra parameter we can set: `shots`. The output of a quantum computing algorithm is probabilistic. Therefore, we must execute the quantum computation multiple times, sampling the outputs to estimate their probabilities. The number of `shots` is the number of executions of the quantum computation. Here we set out `QuantumInstance` to use the least busy backend with 2048 shots.
```python
quantum_instance_real = QuantumInstance(backend_real, shots=2048)
```
The VQE algorithm and QAOA are iterative, meaning that they incorporate a classical-quantum loop which repeats certain computations, _hopefully_ converging to a valid solution. In each iteration, or evaluation, the quantum backend will execute the quantum operations 2048 times. Each shot is quite fast, so we do not have to worry about a significant increase in processing time by using more shots.
----
由於雲端實機免費讓我們練習的 quibits 都很少，我們將問題限制縮小到只需要 4 qubits 來做示範:
組合中角色只包括主教 W 和騎士 M：
1. 主教 W 得 3 分
2. 騎士 M 得 3 分
3. 一位主教 W 搭配一位騎士得 1 分
但是如果派相同角色的人上場會扣分
4. 一位主教 W 和一位主教 W -2 分
5. 一位騎士 M 和一位騎士 M -2 分
用一個二次函數的模型來描述，得分是 $f$:
$f = 3w + 3m + 1wm - 2w^2 - 2m^2$
$0 \le w \le 2$
$0 \le m \le 2$
$w + m \le 4$ (上面 2 個條件就限制了這條，所以這條毋須另外設定)
如果用 Qikit Optimization 套件的工具來描述，則如下:
(程式中變數名稱請取第一個字母來對應題敘)
```python
# Create a small crop-yield example quadratic program
cropyield_small = QuadraticProgram(name="Small Crop-Yield")
# Add two variables, indicating whether we grow 0, 1, or 2 hectares for two different crops
cropyield_small.integer_var(lowerbound=0, upperbound=2, name="Wheat")
cropyield_small.integer_var(lowerbound=0, upperbound=2, name="Maize")
# Add the objective function defining the yield in tonnes
cropyield_small.maximize(
linear={"Wheat": 3, "Maize": 3},
quadratic={("Maize", "Wheat"): 1, ("Maize", "Maize"): -2, ("Wheat", "Wheat"): -2},
)
# This linear constraint is not used as the model never reaches this. This is because the
# sum of the upperbounds on both variables is 4 already. If this constraint is applied, the
# model would require 6 qubits instead of 4.
# cropyield_small.linear_constraint(linear={"Wheat": 1, "Maize": 1}, sense="<=", rhs=4)
print(cropyield_small)
```
\ This file has been generated by DOcplex
\ ENCODING=ISO-8859-1
\Problem name: Small Crop-Yield
Maximize
obj: 3 Wheat + 3 Maize + [ - 4 Wheat^2 + 2 Wheat*Maize - 4 Maize^2 ]/2
Subject To
Bounds
Wheat <= 2
Maize <= 2
Generals
Wheat Maize
End
確定一下確實只需要 4 quibits ...
```python
# Estimate the number of qubits required
ising_operations_small, _ = (
QuadraticProgramToQubo()
.convert(
cropyield_small,
)
.to_ising()
)
print(f"Number of qubits required is {ising_operations_small.num_qubits}")
```
Number of qubits required is 4
接著就可以把問題丟到雲端量子電腦實機解題
```python
# Create our optimizer
optimizer = COBYLA(maxiter=1)
## Get result from real device with VQE
vqe_result_real, vqe_eval_count_real = get_VQE_solution_for(
cropyield_small, quantum_instance_real, optimizer=optimizer
)
```
... (小編補充接續的流程) ...
```python
# Format and print result
print("Solution found using the QAOA method:\n")
print(f"Maximum crop-yield is {vqe_result_real.fval} tons")
print(f"Crops used are: ")
for cropHectares, cropName in zip(vqe_result_real.x, vqe_result_real.variable_names):
print(f"\t{cropHectares} ha of {cropName}")
print(f"\nThe solution was found within {vqe_eval_count_real} evaluations of QAOA.")
```
Solution found using the QAOA method:
Maximum crop-yield is 3.0 tons
Crops used are:
1.0 ha of Wheat
1.0 ha of Maize
The solution was found within 1 evaluations of QAOA.
和傳統解法結果比較一下
```python
# Get classical result
classical_result = get_classical_solution_for(cropyield_small)
# Format and print result
print("Solution found using the classical method:\n")
print(f"Maximum crop-yield is {classical_result.fval} tons")
print(f"Crops used are: ")
_crops = [v.name for v in cropyield_small.variables]
for cropIndex, cropHectares in enumerate(classical_result.x):
print(f"\t{cropHectares} ha of {_crops[cropIndex]}")
```
Solution found using the classical method:
Maximum crop-yield is 3.0 tons
Crops used are:
1.0 ha of Wheat
1.0 ha of Maize
編按：現在我們能試的機器和問題都太小，無法做有意義的效能比較。希望在不久的將來，可以做安兔兔之類的跑分來看看量子計算是不是真的快很多。
## References
[1] A. A. Nel, ‘Crop rotation in the summer rainfall area of South Africa’, South African Journal of Plant and Soil, vol. 22, no. 4, pp. 274–278, Jan. 2005, doi: 10.1080/02571862.2005.10634721.
[2] H. Ritchie and M. Roser, ‘Crop yields’, Our World in Data, 2013, [Online]. Available: https://ourworldindata.org/crop-yields.
[3] G. Brion, ‘Controlling Pests with Plants: The power of intercropping’, UVM Food Feed, Jan. 09, 2014. https://learn.uvm.edu/foodsystemsblog/2014/01/09/controlling-pests-with-plants-the-power-of-intercropping/ (accessed Feb. 15, 2021).
[4] N. O. Ogot, J. O. Pittchar, C. A. O. Midega, and Z. R. Khan, ‘Attributes of push-pull technology in enhancing food and nutrition security’, African Journal of Agriculture and Food Security, vol. 6, pp. 229–242, Mar. 2018.
```python
import qiskit.tools.jupyter
%qiskit_version_table
%qiskit_copyright
```
/Volumes/extra/opt/anaconda3/lib/python3.8/site-packages/qiskit/aqua/__init__.py:86: DeprecationWarning: The package qiskit.aqua is deprecated. It was moved/refactored to qiskit-terra For more information see <https://github.com/Qiskit/qiskit-aqua/blob/main/README.md#migration-guide>
warn_package('aqua', 'qiskit-terra')
<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td><code>qiskit-terra</code></td><td>0.18.1</td></tr><tr><td><code>qiskit-aer</code></td><td>0.8.2</td></tr><tr><td><code>qiskit-ignis</code></td><td>0.6.0</td></tr><tr><td><code>qiskit-ibmq-provider</code></td><td>0.16.0</td></tr><tr><td><code>qiskit-aqua</code></td><td>0.9.4</td></tr><tr><td><code>qiskit</code></td><td>0.29.0</td></tr><tr><td><code>qiskit-nature</code></td><td>0.2.1</td></tr><tr><td><code>qiskit-optimization</code></td><td>0.2.2</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.8.8 (default, Apr 13 2021, 12:59:45)
[Clang 10.0.0 ]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>4</td></tr><tr><td>Memory (Gb)</td><td>8.0</td></tr><tr><td colspan='2'>Sat Sep 18 05:57:06 2021 CST</td></tr></table>
<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>© Copyright IBM 2017, 2021.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>
Sign in with Wallet
Connect another wallet