更新，在自己笔记本上编写的 Latex 公式 在 Github 上渲染显示不对的问题，发现最主要的还是需要添加 反斜杠 转义
Ising模型，伊辛模型
QUBO模型，Quadratic Unconstrained Binary Optimization，二次无约束二进制优化模型
-
自旋变量，其域空间取值是
Ising空间，也就是域空间是$\{-1,1\}$ -
二进制变量，其域空间取值是
布尔空间，也就是域空间是$\{0,1\}$
HOBO，High Order Binary Optimization，高阶二进制优化问题
组合优化问题，组合优化(Combinatorial Optimization, CO)
论文参考： Quantum bridge analytics I: a tutorial on formulating and using QUBO models | SpringerLink
可以将很多组合优化问题重新进行表述为
QUBO模型，而不同类型的约束关系可以用惩罚函数以非常自然的方式体现在“无约束”QUBO公式中。在QUBO公式中，使用惩罚函数可以产生比较精确的模型表示（这句话不是很懂，翻译过来的）。
QUBO模型用于求解NP难问题，该方法用于设计寻找“最优”解决方案的精确求解器基本不可能，除非是非常小的问题实例。该模型使用现代元启发式方法，在有限的时间内找到接近最优解的方案。
最简单的 QUBO 模型可以表述为：
其中
这里的常数方阵
$\mathbf{Q}$ 一般在最优化问题中，在不丧失一般性的情况下，使用 对称矩阵 或者 上/下三角形
对上述进行展开，并改写成标量的形式如下：
由于 QUBO模型表示如下:
其中，
-
$x_i \in \{0,1\}, i \in \mathcal{X}:=\{1, \ldots, X\}$ 是二进制决策变量，并且$\mathcal{E}:=\{(i, j) \mid i, j \in \mathcal{X}, i \neq j\}$ 。 -
$Q_{i j} \in \mathbb{R},(i, j) \in \mathcal{E}$ ,是QUBO目标函数的二次项系数。 -
$c_i \in \mathbb{R}, i \in \mathcal{X}$ ，是QUBO目标函数的一次（线性）项系数。
QUBO问题可以等价的使用Ising模型来表示，通过变量 QUBO模型中决策变量域空间 Ising 模型决策变量域空间 Ising 模型表述如下：
其中，
最优化问题：最小化二元变量的二次目标函数
其中，变量
基于此，将优化模型的最小化目标函数写成矩阵形式：
由于问题规模较小，穷举(
$2^4=16$ 种可能) 求出该问题的最优解为：$(y=-11,x_1=x_4=1,x_2=x_3=0)$
使用python中的
wildqat库（pip install wildqat，这种方法出现错误，有问题的，因为他依赖于python3.6版本的，高版本的 python 库不兼容，安装有问题）import wildqat as wq a = wq.opt() a.qubo = [[-5,2,4,0], [2,-3,1,0], [4,1,-8,5], [0,0,5,-6]] a.sa()这个库有点问题，求的最优解不对，但是可以找到接近最优解。可能是参数设置有问题。
相关代码：
from pyqubo import Binary import neal # 定义哈密顿量 x1, x2, x3, x4 = Binary("x1"), Binary("x2"), Binary("x3"), Binary("x4") H = -5 * x1 - 3 * x2 - 8 * x3 - 6 * x4 + 4 * x1 * x2 + 8 * x1 * x3 + 2 * x2 * x3 + 10 * x3 * x4 # 编译哈密顿量得到一个模型 model = H.compile() # 调用'to_qubo()'获取QUBO系数 # 其中，offset表示下面目标函数中的常数值 # qubo是系数，字典dict类型,{('x2', 'x2'): -3.0,...}表示其系数 qubo, offset = model.to_qubo() print(qubo) print(offset) # 求解qubo模型 # 将Qubo模型输出为BinaryQuadraticModel，BQM来求解 bqm = model.to_bqm() # 定义采样器，这里使用模拟退火采样器 sa = neal.SimulatedAnnealingSampler() # 进行采样，由于该采样器机制，需要设置高一点的采样个数，这样能够确保能够获得最优值 # 对应的设置越高，时间花费越久。 sampleset = sa.sample(bqm, num_reads=10) decoded_samples = model.decode_sampleset(sampleset) # 将采样出来的结果按照目标函数进行排序，求出最佳（小）的采样结果 best_sample = min(decoded_samples, key=lambda x: x.energy) # 输出采样结果 print(f"采样取值：{best_sample.sample}，对应的结果为:{best_sample.energy}")
一般的QUBO模型要求变量为二进制外不包含任何约束。所以如果要实际解决一些有约束条件的问题，需要在目标函数中引入二次惩罚（quadratic penalties）项。
| 经典约束 | 等效惩罚 |
|---|---|
这里注意一下
$P(xy)$ 表示的是标量$P$ 乘以$(x \times y)$ ，也就是$P \times (x \times y)$ ，其他类同。
其中，函数 QUBO中，其约束条件是通过优化器来实现的，因此，惩罚项指定的规则是：对于 最小化问题 的求解，如果是可行解，即满足约束条件，对应的惩罚项等于零；对于不可行解，即不满足约束条件的解，其惩罚项等于一些正的惩罚量。对于
惩罚值太大会阻碍求解过程，因为惩罚项会淹没原始目标函数信息，使得很难区分一个解决方案的质量。另一方面，惩罚值过小会危及寻找可行的解决办法。因此，如何确定惩罚值需要进行思考设计。
对于约束项的设计也有相关的技巧，可以查资料看看。
可以将很多组合优化问题重新进行表述为QUBO模型，并且可以使用模拟退火、量子退火和量子近似优化算法(QAOA)等方法对其求解。都是现代元启发式搜索方法。
很多专有名词和原理细节可以参考：D-Wave System Documentation或者查资料。比如说
纠缠（entanglement）,耦合器（coupler）,本征谱(Eigenspectrum),哈密顿量(Hamiltonian),本征态（eigenstates），初始哈密顿量（Initial Hamiltonian），最终哈密顿量(Final Hamiltonian)，隧穿哈密顿量 (tunneling Hamiltonian)，问题哈密顿量 (problem Hamiltonian)。
量子退火是基于耦合量子位的自然行为来寻找基态（最低能量状态）的过程。量子退火的过程可以使用根据时间变化的哈密顿量
其中，
- 当
$A(0)=1$ 且$B(0)=0$ ,$\mathcal{H}(s)$ 表示初始状态$H_I$ ，由用户自己定义。 - 当
$A(1)=0$ 且$B(1)=0$ ,$\mathcal{H}(s)$ 表示退火后的状态$H_P$ ，这一项也称为问题哈密顿量 (problem Hamiltonian)，是最小能量状态。
那么在设计的时候，初始状态的设计可以根据问题进行简单设计，比如：
其中，
而问题哈密顿量，也就是
其中，
量子退火器首先初始化量子位的叠加态使得
模拟退火算法（Simulated Annealing，SA）是一种通用的全局优化算法，用于在搜索空间中找到最优解。它的基本原理是模拟物理学中的退火过程，通过温度参数控制搜索过程中接受次优解的概率，从而避免陷入局部最优解。 模拟退火算法的基本流程如下：
-
初始化一个解
$x_0$ 和一个初始温度$T_0$ 。 -
在该温度下，进行迭代，每次迭代进行以下操作：
a. 从当前解
$x$ 的邻域中随机选择一个新解$x'$ 。b. 计算
$x'$ 的目标函数值$f(x')$ 和$x$ 的目标函数值$f(x)$ 。c. 根据Metropolis准则判断是否接受新的解，即如果
$f(x')<f(x)$ ，则以概率$1$ 接受$x'$ 作为新的解；否则以概率$e^{-\Delta f/T}$ 接受$x'$ ，其中$\Delta f=f(x')-f(x)$ 。 -
更新温度
$T$ ，并重复步骤 2 直到温度满足停止条件
其中，温度
由上述算法可知，其迭代效果受到初始温度、降温速率和终止温度三个参数的影响。
-
初始温度通常设置为一个较高的值，以保证在搜索的早期能够接受一些劣解，从而有更大的概率跳出局部最优解。
-
终止温度通常设置为一个较小的值，当温度降到终止温度时，搜索停止，此时得到的解就是算法得到的最优解。
-
降温速率通常设置为一个小于1的常数，降温速率越慢，则搜索的时间越长，但是搜索范围也就越广，有更大的概率找到全局最优解。
因此需要合理设置初始温度、降温速率和终止温度三个参数以便能够快速达到最优解。
这里介绍一下dwave-samplers库中的SimulatedAnnealingSampler采样器以及对应的参数。
在该库中，初始温度和终止温度是通过 beta_range ( beta_schedule_type 和 num_sweeps 来设置：对于 beta_schedule_type 参数，有三种方式：
- "linear"：线性，在
$[\beta_0,\beta_1]$ 中线性获取num_sweeps个采样点作为每次更新温度的数值 - "geometric"，几何，在
$[\beta_0,\beta_1]$ 中通过np.geomspace(*beta_range, num=num_betas)获取num_sweeps个采样点作为每次更新温度的数值。 - "custom"：用户自定义
量子近似优化算法可以使用QPanda库，这部分不是很了解，有时间后面补充。
主要解决的问题是：对于QUBO模型中高于二次项的目标函数该如何求解问题
参考论文：[2001.00658] Compressed Quadratization of Higher Order Binary Optimization Problems (arxiv.org)
在不同论文中有不同的叫法：
- Higher Order Binary Optimization (
HOBO)：高阶二进制优化- termwise quadratizationquadratization:按期限平方，这个是将单项式单次转换为二次。也就是
$x_i = x_i^2, \quad where \quad x_i \in \{0,1\}$
高阶二元优化问题的压缩二次化：对于目前而言，存在一种使用 Rosenberg多项式 来 减少高阶优化问题 次数的方法，该方法在布尔空间中通过引入一个额外的变量来减少一个项的次数，这种方法最终会导致最后的布尔空间矩阵很稀疏，并且问题规模会变大。然后在论文中，提出了一种直接在原有布尔空间中运行的降阶方法，而不通过增加额外变量的方式。
论文中有个没看懂：
感觉这里的
$x$ 打错了，应该是$z$ ，然后想表达的意思应该是，对于$z$ 属于布尔空间的取值，可以通过增加变量$y$ ，使得当$y$ 取某个 特定值 的时候$h(z,y_{取特定值})$ 等价于$f(z)$ ，由于$h$ 函数不确定，只是说能找到该函数满足等价这个条件，那么再做一个限定，也就是当$y$ 取该 特定值 的时候，正好是$h(z,y_{取特定值})$ 取到最小值，那么根据贪心思想，很容易理解，到时候 最小化
$f(z)$ 等价于 最小化$h(z,y)$
基于此，就有等价：
这个比较难理解，但是，想到 Linear and quadratic reformulations of nonlinear optimization problems in binary variables中证明了对于最高次数
这种代替方式有两个问题，具体可以看看论文，我们没有使用变量代替降次，所以对于这个
HOBO没有仔细研究。
| 名词 | 数学表达 |
|---|---|
| 贷款金额（Loan amount） |
|
| 利息收入率（Interest income rate） |
|
| 通过率矩阵， |
其中， |
| 坏账率矩阵， |
其中， |
| 总通过率 | |
| 总坏账率 | |
| 集合空间，具体在文章中说明 | |
| 约束项 | |
| 最终收入 |
对于问题一，我们定义决策变量
那么对于决策变量
于是决策变量
考虑整体，也就是最终收入
对于最终收入，我们考虑最小化目标函数
考虑约束条件：由于在整体中只有一个评分卡的阈值能被选取，因此在决策变量
这里的
$\mathcal{X}$ 空间定义为：$\mathcal{X}:= { (i, j) \mid i \in {1,2,...,100 }, j \in { 1,2,...,10 } }$
综上所述：
通过添加约束项 QUBO）:
对于问题一的 枚举暴力 代码实现：
import time
import pandas as pd
if __name__ == '__main__':
# 定义一些变量
L = 1000000 # 贷款资金
I = 0.08 # 利息收入率
df = pd.read_csv('./data/附件1：data_100.csv', header=0)
# 定义最大的数
max_value = -1
# 定义最大数字的索引下标
max_i = max_j = -1
start = time.time()
for i in range(1, 101):
for j in range(0, 10):
# 总通过率
P = df[f"t_{i}"].iloc[j]
# 总坏账率
Q = df[f"h_{i}"].iloc[j]
# 此时的最终收入
temp_value = L * I * P * (1 - Q) - L * P * Q
# print(f"选卡[{i},{j + 1}]，对应最终收入{temp_value},最大收入:{max_value}")
if temp_value > max_value:
max_value = temp_value
max_i = i
max_j = j
end = time.time()
print(f"暴力时间：{end - start} s")
print(f"最大值：{max_value}")
print(f"对应的选取卡1：第{max_i}张{max_j + 1}个阈值")使用 QUBO 模型+模拟退火算法实现：
import time
from pyqubo import Array, Constraint
from dwave.samplers import SimulatedAnnealingSampler
import numpy as np
from pyqubo import Placeholder
# 在字典中找到value的key
def get_key(dict, value):
return [k for k, v in dict.items() if v == value]
# 主函数,问题一的主要求解函数
def main():
# 定义一些变量
L = 1000000 # 贷款资金
I = 0.08 # 利息收入率
# 读取数据
data = np.genfromtxt('./data/附件1：data_100.csv', delimiter=',', skip_header=1)
# 获取T矩阵和H矩阵
T = data[:, ::2]
H = data[:, 1::2]
# 定义二进制决策变量,10 * 100
x = Array.create('x', shape=(10, 100), vartype='BINARY')
# 定义惩罚项M
M = Placeholder('M')
M = 100000
# 定义哈密顿量,也就是对应的目标函数
# 注意这里不是总通过率和总坏账率，是中间变量，将最终决策目标的连加符号放到里面整出来的
P = np.sum(np.multiply(x, T))
Q = np.sum(np.multiply(x, H))
H = - (L * I * P * (1 - Q) - L * P * Q) + M * Constraint((np.sum(x) - 1) ** 2, label='sum(x_i_j) = 1')
# 编译哈密顿量得到一个模型
model = H.compile()
# 将Qubo模型输出为BinaryQuadraticModel，BQM来求解
bqm = model.to_bqm()
# 记录开始退火时间
start = time.time()
# 模拟退火
sa = SimulatedAnnealingSampler()
sampleset = sa.sample(bqm, seed=666, beta_range=[10e-10, 50], num_sweeps=10000, beta_schedule_type='geometric',num_reads=10)
# 对数据进行筛选，对选取的数据进行选取最优的
decoded_samples = model.decode_sampleset(sampleset) # 将上述采样最好的num_reads组数据变为模型可读的样本数据
best_sample = min(decoded_samples, key=lambda x: x.energy) # 将能量值最低的样本统计出来，表示BQM的最优解
end = time.time()
print(f"退火时间花费：{end - start} s")
# 统计决策变量为1的所有数据
data_1_list = get_key(best_sample.sample, 1)
print(f"对应的取1的决策变量有：{data_1_list}(注意是索引，从0开始),对应的能量为(最大最终收入):{- best_sample.energy}")
if __name__ == '__main__':
main()对于问题一，虽然暴力穷举比
QUBO模型+模拟退火算法运行时间快，但是考虑算法时间复杂度：卡有$m$ 张，阈值选取有$n$ 个
- 暴力穷举算法，由于有循环，因此是
$O(mn)$ QUBO模型+模拟退火算法，属于是现代元启发式优化算法，在数据较多时，尤其是NP问题上，规模越大，其运行速度相对暴力要快的多。
对于问题二，我们考虑了两种QUBO模型的理解：
- 定义决策变量
$x_{i,j,k}$ 表示是否选：第1张信用评分卡选第$i$ 个阈值，第2张信用评分卡选第$j$ 个阈值，第3张信用评分卡选第$k$ 个阈值。 - 主要来源于暴力穷举过程的一种发现，固定前两张信用评分卡的阈值选择，定义第三张信用评分卡第
$i$ 个阈值是否选择作为决策变量$x_i$ ，决策出第三张信用评分卡的最佳阈值后，同理决策出其他两张信用卡的最佳阈值。多迭代几轮，这种方法决策出来的阈值在该问题中就是最佳的阈值组合。
具体看下面内容：
定义决策变量
那么对于决策变量
于是决策变量
考虑整体，也就是最终收入
同样考虑最小化目标函数
这里的
$\mathcal{X}$ 空间定义为：$\mathcal{X}:={(i,j,k) \mid i,j,k \in {1,2,...,10} }$
综上所述：
通过添加约束项 QUBO）:
在代码实现过程中，由于需要计算总通过率
$T_{1,i} \ * \ T_{2,j} \ * \ T_{3,k}$ 以及总坏账率$\frac{1}{3} * (H_{1,i} + H_{2,j} + H_{3,k})$ 。如果使用for循环效率比较慢，这里使用了numpy库中的广播机制来获得相关计算矩阵。
代码实现如下：
import time
from pyqubo import Array, Constraint
from dwave.samplers import SimulatedAnnealingSampler
import numpy as np
from pyqubo import Placeholder
# 在字典中找到value的key
def get_key(dict, value):
return [k for k, v in dict.items() if v == value]
def main():
# 定义一些变量
L = 1000000 # 贷款资金
I = 0.08 # 利息收入率
# 表示选取的卡号
card1, card2, card3 = 1, 2, 3
# 读取数据
data = np.genfromtxt('../data/附件1：data_100.csv', delimiter=',', skip_header=1)
# 获取T矩阵和H矩阵
T = data[:, ::2]
H = data[:, 1::2]
# 定义二进制决策变量,10 * 100
x = Array.create('x', shape=(10, 10, 10), vartype='BINARY')
# 定义惩罚项M
M = Placeholder('M')
M = 30000
# 定义哈密顿量,也就是对应的目标函数
# 注意这里不是总通过率和总坏账率，是中间变量，将最终决策目标的连加符号放到里面整出来的
T1 = T[:, card1 - 1]
T2 = T[:, card2 - 1]
T3 = T[:, card3 - 1]
H1 = H[:, card1 - 1]
H2 = H[:, card2 - 1]
H3 = H[:, card3 - 1]
# 计算三种组合后的总通过率和总坏账率,通过广播机制来实现
T = T1[:, None, None] * T2[None, :, None] * T3[None, None, :]
H = (H1[:, None, None] + H2[None, :, None] + H3[None, None, :]) / 3
P = np.sum(np.multiply(x, T))
Q = np.sum(np.multiply(x, H))
H = - (L * I * P * (1 - Q) - L * P * Q) + M * Constraint((np.sum(x) - 1) ** 2, label='sum(x_i_j) = 1')
# 编译哈密顿量得到一个模型
model = H.compile()
# 将Qubo模型输出为BinaryQuadraticModel，BQM来求解
bqm = model.to_bqm()
print("开始模拟退火")
start = time.time()
sa = SimulatedAnnealingSampler()
sampleset = sa.sample(bqm, seed=888, beta_range=[10e-12, 60], beta_schedule_type='geometric', num_reads=50)
# 对数据进行筛选
decoded_samples = model.decode_sampleset(sampleset) # 将上述采样最好的num_reads组数据变为模型可读的样本数据
best_sample = min(decoded_samples, key=lambda x: x.energy) # 将能量值最低的样本统计出来，表示BQM的最优解
# print(f"验证约束条件M：{best_sample.constraints()}")
end = time.time()
print(f"退火时间：{end - start} s")
# print(best_sample.sample)
# 统计决策变量为1的所有数据
data_1_list = get_key(best_sample.sample, 1)
print(f"对应的取1的决策变量有：{data_1_list}(表示[第一张卡的阈值][第二张卡的阈值][第三张卡的阈值]，注意是索引，从0开始),对应的能量为{-best_sample.energy}")
if __name__ == '__main__':
main()对应的枚举暴力代码如下：
import pandas as pd
if __name__ == '__main__':
# 定义一些变量
L = 1000000 # 贷款资金
I = 0.08 # 利息收入率
# 选取第几列的信用卡
card1 = 1
card2 = 2
card3 = 3
df = pd.read_csv('../data/附件1：data_100.csv', header=0)
# 定义最大的数
max_value = -1
# 定义最大数字的索引下标
max_i = max_j = max_k = -1
for i in range(0, 10):
for j in range(0, 10):
for k in range(0, 10):
# 总通过率
P = df[f"t_{card1}"].iloc[i] * df[f"t_{card2}"].iloc[j] * df[f"t_{card3}"].iloc[k]
# 总坏账率
Q = (df[f"h_{card1}"].iloc[i] + df[f"h_{card2}"].iloc[j] + df[f"h_{card3}"].iloc[k]) / 3
# 此时的最终收入
temp_value = L * I * A * (1 - B) - L * A * B
# print(f"选卡1阈值：[{i + 1}]，选卡1阈值：[{j + 1}]，选卡1阈值：[{k + 1}]，对应最终收入{temp_value},之前的最大收入:{max_value}")
if temp_value > max_value:
max_value = temp_value
max_i = i
max_j = j
max_k = k
print(f"第{i + 1}个对应的最大值：{max_value},三个阈值选取为:{max_i + 1},{max_j + 1},{max_k + 1},其中，阈值编号[1-10]")
print(f"最大值：{max_value},三个阈值选取为:{max_i + 1},{max_j + 1},{max_k + 1},其中，阈值编号[1-10]")在运行过程中，发现一个对第三问很有帮助的规律：可以把上面的决策过程（for循环）打印出来看一下：
对于第一个循环，其决策过程是：[1,1,2]-->[2,1,2]-->[3,1,2]-->[3,1,2]-->[3,1,2]-->[6,1,2]-->[7,1,2]-->[8,1,2]-->[8,1,2]-->[8,1,2]。这里打印出来结果比较少，可能看的不是很明显规律。但是容易发现：
对于三张评分信用卡的决策过程，如果固定第一张和第二张评分信用卡的阈值选取，而考虑第三张评分信用卡的阈值选取，此时选择最优的信用评分卡阈值。然后固定第一张和第三张评分信用卡的阈值选取，而考虑第二张评分信用卡的阈值选取，此时也选择最优的信用评分卡阈值。然后固定第二张和第三张评分信用卡的阈值选取，而考虑第一张评分信用卡的阈值选取，此时也选择最优的信用评分卡阈值。基于此不停的迭代，那么最终获得的阈值选取就是最优的阈值组合。（可以考虑证明一下）。
基于该过程，我们就可以将原本
最开始我们考虑的是：设计决策变量
在不考虑约束条件的情况下，整体的通过率
对应的最终收入
可以发现，如果使用
如果使用
那么基于这个理论就可以将原本的目标函数转换为最高二次项的目标函数。不过这个也有点问题，就是由于决策变量
以上所有考虑都存在一个问题，就是怎样解决整体的通过率
那么对于决策变量
那么带入最终收入公式，有：
这种方法如果不考虑约束条件，那么量子个数有
考虑使用循环暴力遍历其找到最佳的组合。这个需要很长的时间来进行，那么考虑多线程+分治的思想对其暴力，代码如下：
import time
import pandas as pd
from threading import Thread
# 定义一些变量
L = 1000000 # 贷款资金
I = 0.08 # 利息收入率
df = pd.read_csv('../data/附件1：data_100.csv', header=0)
# 全局变量，用于存储所有线程中的最大值
max_value = -1
# 定义最大数字的索引下标
max_i = max_j = max_k = max_m = max_n = max_l = -1
def computer_result(begin, end):
# 使用全局变量
global max_value, max_i, max_j, max_k, max_m, max_n, max_l
for i in range(begin, end):
for j in range(i + 1, 101):
for k in range(j + 1, 101):
for m in range(0, 10):
for n in range(0, 10):
for l in range(0, 10):
# 总通过率
A = df[f"t_{i}"].iloc[m] * df[f"t_{j}"].iloc[n] * df[f"t_{k}"].iloc[l]
# 总坏账率
B = (df[f"h_{i}"].iloc[m] + df[f"h_{j}"].iloc[n] + df[f"h_{k}"].iloc[l]) / 3
# 此时的最终收入
temp_value = L * I * A * (1 - B) - L * A * B
if temp_value > max_value:
max_value, max_i, max_j, max_k, max_m, max_n, max_l = temp_value, i, j, k, m, n, l
if __name__ == '__main__':
start = time.time()
print(start)
# 创建线程并启动它们
threads = [Thread(target=computer_result, args=(i, i + 1)) for i in range(1, 101)]
for thread in threads:
thread.start()
# 等待所有线程执行完毕
for thread in threads:
thread.join()
end = time.time()
print(f"时间花费：{end - start} s")
print(f"最大值：{max_value},对应的选取卡1：第{max_i}张{max_m + 1}个阈值,"
f"对应的选取卡2：第{max_j}张{max_n + 1}个阈值,"
f"对应的选取卡3：第{max_k}张{max_l + 1}个阈值,其中，卡取值[1-100],阈值取值[1-10]")感觉代码还是可以优化的，不过不太想花时间在这个上面。
主要是问题二中找到的贪心规律，所以这个效率比较高。那么类似的使用如下迭代算法：
那么最主要的问题就是固定两个评分信用卡的阈值后，如何找到最优的信用评分卡的阈值？要么暴力，要么使用问题一的QUBO。
import pandas as pd
if __name__ == '__main__':
# 定义一些变量
L = 1000000 # 贷款资金
I = 0.08 # 利息收入率
# 迭代次数
Iter = 10
# 选取第几列的信用卡,初始化最开始选择的数据
card1, threshold1, card2, threshold2, card3, threshold3 = 1, 0, 2, 0, 3, 0
# 定义最大的数
max_value = -1
# 读取数据
df = pd.read_csv('../data/附件1：data_100.csv', header=0)
for iter in range(Iter): # 迭代次数，其实这里还可以根据几次迭代算出的总值来提前结束迭代
# 设置更新迭代,不使用临时变量
card1, threshold1, card2, threshold2, card3, threshold3 = card2, threshold2, card3, threshold3, card1, threshold1
# 从固定的信用卡1和信用卡2外数据中选一个，构建可以选择的卡的集合
choose_list = {i for i in range(1, 101) if i not in (card1, card2)}
# 获取第一张卡和第二张卡的通过率和坏账率，否则在循环里重复获取了
t1 = df[f"t_{card1}"].iloc[threshold1]
t2 = df[f"t_{card2}"].iloc[threshold2]
h1 = df[f"h_{card1}"].iloc[threshold1]
h2 = df[f"h_{card2}"].iloc[threshold2]
# 固定第一张卡,和第二张卡，然后找到最佳的第三张卡
for c3 in choose_list:
for th3 in range(0, 10):
# 总通过率
A = t1 * t2 * df[f"t_{c3}"].iloc[th3]
# 总坏账率
B = (h1 + h2 + df[f"h_{c3}"].iloc[th3]) / 3
# 此时的最终收入
temp_value = L * I * A * (1 - B) - L * A * B
if temp_value > max_value:
max_value = temp_value
card3 = c3
threshold3 = th3
print(f"最大值：{max_value},卡选择:{card1}_{threshold1 + 1},{card2}_{threshold2 + 1},{card3}_{threshold3 + 1}")代码：
import pandas as pd
import time
import re
from pyqubo import Array, Constraint
from dwave.samplers import SimulatedAnnealingSampler
import numpy as np
from pyqubo import Placeholder
# 在字典中找到value的key
def get_key(dict, value):
return [k for k, v in dict.items() if v == value]
def main():
# 定义一些变量
L = 1000000 # 贷款资金
I = 0.08 # 利息收入率
# 迭代次数
Iter = 10
# 选取第几列的信用卡,初始化最开始选择的数据
card1, threshold1, card2, threshold2, card3, threshold3 = 1, 0, 2, 0, 3, 0
# 定义最大的数
max_value = -1
# 读取数据
df = pd.read_csv('../data/附件1：data_100.csv', header=0)
for iter in range(Iter): # 迭代次数，其实这里还可以根据几次迭代算出的总值来提前结束迭代
# 设置更新迭代,不使用临时变量
card1, threshold1, card2, threshold2, card3, threshold3 = card2, threshold2, card3, threshold3, card1, threshold1
# 获取第一张卡和第二张卡的通过率和坏账率
t1 = df[f"t_{card1}"].iloc[threshold1]
t2 = df[f"t_{card2}"].iloc[threshold2]
h1 = df[f"h_{card1}"].iloc[threshold1]
h2 = df[f"h_{card2}"].iloc[threshold2]
# 获取第三张卡能够选取的通过率与坏账率的矩阵
choose_index = [i for i in range(1, 101) if i not in (card1, card2)]
card3_data = df.drop([f"t_{card1}", f"h_{card1}", f"t_{card2}", f"h_{card2}"], axis=1).values
# 获取T矩阵和H矩阵
T = card3_data[:, ::2]
H = card3_data[:, 1::2]
# 定义二进制决策变量,10 * 98
x = Array.create('x', shape=(10, 98), vartype='BINARY')
# 定义惩罚项M
M = Placeholder('M')
M = 50000
# 定义哈密顿量,也就是对应的目标函数
P = t1 * t2 * np.sum(np.multiply(x, T))
Q = (h1 + h2 + np.sum(np.multiply(x, H))) / 3
H = - (L * I * P * (1 - Q) - L * P * Q) + M * Constraint((np.sum(x) - 1) ** 2, label='sum(x_i_j) = 1')
# 编译哈密顿量得到一个模型
model = H.compile()
bqm = model.to_bqm()
# 记录开始退火时间
start = time.time()
#
sa = SimulatedAnnealingSampler()
sampleset = sa.sample(bqm, seed=666, beta_range=[10e-10, 50], num_sweeps=999, beta_schedule_type='geometric',
num_reads=50)
# 对数据进行筛选，对选取的数据进行选取最优的
decoded_samples = model.decode_sampleset(sampleset) # 将上述采样最好的num_reads组数据变为模型可读的样本数据
best_sample = min(decoded_samples, key=lambda x: x.energy) # 将能量值最低的样本统计出来，表示BQM的最优解
end = time.time()
# print(f"退火时间花费：{end - start} s")
# 统计决策变量为1的所有数据,并对第一个数据进行拆解，获得对应的下标
data_1_list = get_key(best_sample.sample, 1)
index_i_j = re.findall(r'\d+', data_1_list[0])
index_i = int(index_i_j[0])
if max_value < - best_sample.energy:
card3 = choose_index[int(index_i_j[1])]
threshold3 = index_i
max_value = - best_sample.energy
print(
f"第{iter + 1}次迭代，最大值：{max_value},卡选择:{card1}_{threshold1 + 1},{card2}_{threshold2 + 1},{card3}_{threshold3 + 1}")
print(f"最大值：{max_value},卡选择:{card1}_{threshold1 + 1},{card2}_{threshold2 + 1},{card3}_{threshold3 + 1}")
if __name__ == '__main__':
start = time.time()
main()
end = time.time()
print(f"时间花费： {end - start}s")