一些名词

Ising模型，伊辛模型

QUBO模型，Quadratic Unconstrained Binary Optimization，二次无约束二进制优化模型

1.自旋变量，其域空间取值是Ising空间，也就是域空间是{ − 1 , 1 } 

2.二进制变量，其域空间取值是布尔空间，也就是域空间是{ 0 , 1 } 

HOBO，High Order Binary Optimization，高阶二进制优化问题

组合优化问题，组合优化(Combinatorial Optimization, CO)

二次无约束二进制优化(QUBO)模型

论文参考： Quantum bridge analytics I: a tutorial on formulating and using QUBO models | SpringerLink

可以将很多组合优化问题重新进行表述为QUBO模型，而不同类型的约束关系可以用惩罚函数以非常自然的方式体现在“无约束”QUBO公式中。在QUBO公式中，使用惩罚函数可以产生比较精确的模型表示（这句话不是很懂，翻译过来的）。

QUBO模型用于求解NP难问题，该方法用于设计寻找“最优”解决方案的精确求解器基本不可能，除非是非常小的问题实例。该模型使用现代元启发式方法，在有限的时间内找到接近最优解的方案。

QUBO基本模型

最简单的QUBO模型可以表述为：

其中 x是二进制决策变量的向量，Q 是常数的方阵

这里的常数方阵Q 一般在最优化问题中，在不丧失一般性的情况下，使用对称矩阵或者上/下三角形

对上述进行展开，并改写成标量的形式如下：

由于x _i ∈ { 0 , 1 } 属于二进制变量，存在x _i =，因此对于更具一般性的QUBO模型表示如下:

其中，

QUBO问题可以等价的使用Ising模型来表示，通过变量y_i​=1−2x_i​进行转换，将原本QUBO模型中决策变量域空间{ 0 , 1 }映射到Ising模型决策变量域空间{ − 1 , + 1 }。转变后的Ising模型表述如下：

其中，y_i​∈{−1,1},i∈X

论文中一个小小的案例

最优化问题：最小化二元变量的二次目标函数 y：

其中，变量x_i​∈{0,1}。它有一个线性部分（−5x₁−3x₂−8x₃−6x₄ ）和一个二次部分（4x₁ x₂ +8x₁ x₃ +2x₂ x₃ +10x₃ x₄），由于二进制变量x _j 满足x _j =，所以线性部分可以写成：

基于此，将优化模型的最小化目标函数写成矩阵形式：

由于问题规模较小，穷举(2⁴=16 种可能) 求出该问题的最优解为：y = − 11 , x ₁ = x ₄ = 1 , x ₂ = x ₃ = 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()

这个库有点问题，求的最优解不对，但是可以找到接近最优解。可能是参数设置有问题。

 

使用PyQUBO库：[2103.01708] PyQUBO: Python Library for Mapping Combinatorial Optimization Problems to QUBO Form (arxiv.org)

相关代码：

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模型

一般的QUBO模型要求变量为二进制外不包含任何约束。所以如果要实际解决一些有约束条件的问题，需要在目标函数中引入二次惩罚（quadratic penalties）项。

 这里注意一下P (xy) 表示的是标量P  乘以( x ∗ y )，也就是P ∗ ( x ∗ y )，其他类同。

其中，函数 P计算出来的是正的标量惩罚值，且这个值计算的必须得足够大。因为，在QUBO中，其约束条件是通过优化器来实现的，因此，惩罚项指定的规则是：对于最小化问题的求解，如果是可行解，即满足约束条件，对应的惩罚项等于零；对于不可行解，即不满足约束条件的解，其惩罚项等于一些正的惩罚量。对于P 的取值是否合理可以通过求解模型的解，去计算添加的惩罚函数，如果是0，那就说明，P值效果可以。

惩罚值太大会阻碍求解过程，因为惩罚项会淹没原始目标函数信息，使得很难区分一个解决方案的质量。另一方面，惩罚值过小会危及寻找可行的解决办法。因此，如何确定惩罚值需要进行思考设计。

对于约束项的设计也有相关的技巧，可以查资料看看。

求解

可以将很多组合优化问题重新进行表述为QUBO模型，并且可以使用模拟退火、量子退火和量子近似优化算法(QAOA)等方法对其求解。都是现代元启发式搜索方法。

量子退火

很多专有名词和原理细节可以参考：D-Wave System Documentation或者查资料。比如说纠缠（entanglement）,耦合器（coupler）,本征谱(Eigenspectrum),哈密顿量(Hamiltonian),本征态（eigenstates），初始哈密顿量（Initial Hamiltonian），最终哈密顿量(Final Hamiltonian)，隧穿哈密顿量 (tunneling Hamiltonian)，问题哈密顿量 (problem Hamiltonian)。

量子退火是基于耦合量子位的自然行为来寻找基态（最低能量状态）的过程。量子退火的过程可以使用根据时间变化的哈密顿量来表示：

其中， A ( s ) 和B ( s )是退火路径函数，根据归一化退火时间s = t / t _a ，其中t _a 表示总的退火时间。这样设计可以保证：

那么在设计的时候，初始状态的设计可以根据问题进行简单设计，比如：

其中，表示第 i个量子位的P a u l i − x 算子（泡利− x 算子）

而问题哈密顿量，也就是H _P设计可通过最优化的目标函数给出：

其中，表示第 i个量子位的P a u l i − z 算子（泡利− z算子）。

量子退火器首先初始化量子位的叠加态使得，然后在退火时间内操作量子位耦合，使朝向 H _P 发展，退火结束后，就处于 problem Hamiltonian 的本征态（eigenstate）。根据量子计算的绝热定理，如果退火时间足够长，该时变哈密顿量将始终保持基态，即最优化函数（的解）

模拟退火

模拟退火算法（Simulated Annealing，SA）是一种通用的全局优化算法，用于在搜索空间中找到最优解。它的基本原理是模拟物理学中的退火过程，通过温度参数控制搜索过程中接受次优解的概率，从而避免陷入局部最优解。 模拟退火算法的基本流程如下：

1.初始化一个解 x ₀和一个初始温度 T ₀ 。

2.在该温度下，进行迭代，每次迭代进行以下操作：

a.从当前解 x的邻域中随机选择一个新解 x ′。

b.计算 x ′的目标函数值 f ( x ′ ) 和 x的目标函数值 f ( x )。

c.根据Metropolis准则判断是否接受新的解，即如果 f ( x ′ ) < f ( x )，则以概率1接受 x ′  作为新的解；否则以概率 e − ^{Δ f / T} 接受 x ′ ，其中 Δ f = f ( x ′ ) − f (x) 。

3.更新温度 T，并重复步骤 2 直到温度满足停止条件

其中，温度 T 是一个随着时间逐渐降低的参数，控制着算法接受次优解的概率。

由上述算法可知，其迭代效果受到初始温度、降温速率和终止温度三个参数的影响。

1.初始温度通常设置为一个较高的值，以保证在搜索的早期能够接受一些劣解，从而有更大的概率跳出局部最优解。

2.终止温度通常设置为一个较小的值，当温度降到终止温度时，搜索停止，此时得到的解就是算法得到的最优解。

3.降温速率通常设置为一个小于1的常数，降温速率越慢，则搜索的时间越长，但是搜索范围也就越广，有更大的概率找到全局最优解。

因此需要合理设置初始温度、降温速率和终止温度三个参数以便能够快速达到最优解。

这里介绍一下dwave-samplers库中的SimulatedAnnealingSampler采样器以及对应的参数。

在该库中，初始温度和终止温度是通过转换为β ₁和β ₀ （其中，K是波尔兹常量），通过传入参数 beta_range = [ β 0 , β 1 ] 传入。而降温速率一般通过beta_schedule_type和num_sweeps来设置：对于beta_schedule_type参数，有三种方式：

1. “linear”：线性，在[ β 0 , β 1 ]中线性获取num_sweeps个采样点作为每次更新温度的数值
2. “geometric”，几何，在[ β 0 , β 1 ]中通过np.geomspace(*beta_range, num=num_betas)获取num_sweeps个采样点作为每次更新温度的数值。
3. “custom”：用户自定义

量子近似优化算法(QAOA)

量子近似优化算法可以使用QPanda库，这部分不是很了解，有时间后面补充。

高阶问题（HOBO）

主要解决的问题是：对于QUBO模型中高于二次项的目标函数该如何求解问题

参考论文：[2001.00658] Compressed Quadratization of Higher Order Binary Optimization Problems (arxiv.org)

在不同论文中有不同的叫法：

Higher Order Binary Optimization (HOBO)：高阶二进制优化
termwise quadratizationquadratization:按期限平方，这个是将单项式单次转换为二次。也就是x _i =   , where x _i ∈ 0 , 1 

高阶二元优化问题的压缩二次化：对于目前而言，存在一种使用Rosenberg多项式来减少高阶优化问题次数的方法，该方法在布尔空间中通过引入一个额外的变量来减少一个项的次数，这种方法最终会导致最后的布尔空间矩阵很稀疏，并且问题规模会变大。然后在论文中，提出了一种直接在原有布尔空间中运行的降阶方法，而不通过增加额外变量的方式。

论文中有个没看懂：

感觉这里的x打错了，应该是z，然后想表达的意思应该是，对于z属于布尔空间的取值，可以通过增加变量 y，使得当 y取某个特定值的时候 h ( _{z , y 取特定值} )等价于f ( z )，由于 h函数不确定，只是说能找到该函数满足等价这个条件，那么再做一个限定，也就是当y 取该特定值的时候，正好是h ( _{z , y 取特定值} )取到最小值

那么根据贪心思想，很容易理解，到时候最小化f ( z )等价于最小化 h ( z , y ) 

基于此，就有等价：

这个比较难理解，但是，想到x _i ∈ { 0 , 1 }，并且，这里y是很多个的组合，并且y _i ∈ { 0 , 1 } 就比较容易穷举，在论文Linear and quadratic reformulations of nonlinear optimization problems in binary variables中证明了对于最高次数d ，需要的变量y个数是⌈ log ⁡ d ⌉ − 1。（感觉公式表达有问题）

这种代替方式有两个问题，具体可以看看论文，我们没有使用变量代替降次，所以对于这个HOBO没有仔细研究。

2023年mathorcup的A题

数学定义

问题一

对于问题一，我们定义决策变量 x _{i , j​}用于表示：

那么对于决策变量 x _{i , j}对应的通过率 P _{i , j} 与坏损率 Q _{i , j}有：

于是决策变量 x _{i , j} 对应的最终收入 H _{i , j} 有

考虑整体，也就是最终收入 H有

对于最终收入，我们考虑最小化目标函数 − H ,即

考虑约束条件：由于在整体中只有一个评分卡的阈值能被选取，因此在决策变量x _{i , j} 中只有一个是1，即：

这里的 X空间定义为：X : = { ( i , j ) | i ∈ { 1 , 2 , . . . , 100 } , j ∈ { 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 ，决策出第三张信用评分卡的最佳阈值后，同理决策出其他两张信用卡的最佳阈值。多迭代几轮，这种方法决策出来的阈值在该问题中就是最佳的阈值组合。

具体看下面内容：

QUBO

定义决策变量 x _{i , j , k}用于表示:

x _{i , j , k=}

那么对于决策变量 x _{i , j , k},对应的通过率 P _{i , j , k}与坏损率 Q _{i , j , k}有：

于是决策变量 x _{i , j , k} 对应的最终收入 H _{i , j , k} 有

考虑整体，也就是最终收入 H有

同样考虑最小化目标函数 − H ,此时决策变量的约束条件为：(只有一种方案选取)

这里的 X空间定义为：X : = { ( i , j , k ) |  i , j , k ∈ { 1 , 2 , . . . , 10 } } 

综上所述：

通过添加约束项 λ 将上述二次有约束二进制优化模型改写为二次无约束二进制优化模型（QUBO）:

在代码实现过程中，由于需要计算总通过率T _{1 , i}  ∗  T _{2 , j}   ∗  T _{3 , k}以及总坏账率 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 * df[f"t_{card2}"].iloc[j] * df[f"t_{card3}"].iloc[k]

                # 总坏账率
                Q = (df[f"h_{card1}"].iloc + df[f"h_{card2}"].iloc[j] + df[f"h_{card3}"].iloc[k]) / 3

                # 此时的最终收入
                temp_value = L * I * P * (1 - Q) - L * P * Q

                # 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]")

贪心+ QUBO

在运行过程中，发现一个对第三问很有帮助的规律：可以把上面的决策过程（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]。这里打印出来结果比较少，可能看的不是很明显规律。但是容易发现：

对于三张评分信用卡的决策过程，如果固定第一张和第二张评分信用卡的阈值选取，而考虑第三张评分信用卡的阈值选取，此时选择最优的信用评分卡阈值。然后固定第一张和第三张评分信用卡的阈值选取，而考虑第二张评分信用卡的阈值选取，此时也选择最优的信用评分卡阈值。然后固定第二张和第三张评分信用卡的阈值选取，而考虑第一张评分信用卡的阈值选取，此时也选择最优的信用评分卡阈值。基于此不停的迭代，那么最终获得的阈值选取就是最优的阈值组合。（可以考虑证明一下）。

基于该过程，我们就可以将原本 10 × 10 × 10 个决策变量的求解转换为 迭代次数 × 10 个决策变量的问题。这对于问题三而说是非常重要的一种发现。对于问题三，原本是 1000 × 990 × 980个决策变量问题转换为了 迭代次数 × 980 个决策变量规模的问题。

问题三

最开始我们考虑的是：设计决策变量x _{i , j} 用于表示：

在不考虑约束条件的情况下，整体的通过率P和坏损率Q 满足：

对应的最终收入H有：

可以发现，如果使用这个公式，其目标函数中，决策变量的最高次项系数为1001次。转换为了 HOBO模型，需要降次，感觉很复杂。不考虑。

如果使用这个公式，然后考虑使用多项式近似(Polynomial approximation)的思想，将指数项使用多项式展开，比如泰勒级数展开式：,保留一阶项，即：。也就可以将指数项转换为一个多项式函数的线性项：

那么基于这个理论就可以将原本的目标函数转换为最高二次项的目标函数。不过这个也有点问题，就是由于决策变量x _{i , j} ∈ { 0 , 1 } 是离散变量，那么使用泰勒级展开式未必等价， 也不做考虑。

以上所有考虑都存在一个问题，就是怎样解决整体的通过率P表述。如果能够将P 转换为一次项的表达，那么对于整个求解问题就很简单了。那么就重新定义决策变量x _{i , j , k , l , m , n}用于决策是否选择：第 i张信用评分卡选第 l 个阈值，第 j张信用评分卡选第 m 个阈值，第 k 张信用评分卡选第 n个阈值。即：

那么对于决策变量 x _{i , j , k , l , m , n} ，对应的通过率 P _{i , j , k , l , m , n} 与坏损率 Q _{i , j , k , l , m , n} 有：

那么对于决策变量 x _{i , j , k , l , m , n} 对应的通过率 P _{i , j , k , l , m , n} 与坏损率 Q _{i , j , k , l , m , n} 有：

那么代入最终收入公式，有：

=

这种方法如果不考虑约束条件，那么量子个数有 100 × 100 × 100 × 10 × 10 × 10 ，内存有点不太够，运行不出来这个。

考虑暴力

考虑使用循环暴力遍历其找到最佳的组合。这个需要很长的时间来进行，那么考虑多线程+分治的思想对其暴力，代码如下：

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):
                            # 总通过率
                            P = df[f"t_{i}"].iloc[m] * df[f"t_{j}"].iloc[n] * df[f"t_{k}"].iloc[l]

                            # 总坏账率
                            Q = (df[f"h_{i}"].iloc[m] + df[f"h_{j}"].iloc[n] + df[f"h_{k}"].iloc[l]) / 3

                            # 此时的最终收入
                            temp_value = L * I * P * (1 - Q) - L * P * Q

                            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):
                # 总通过率
                P = t1 * t2 * df[f"t_{c3}"].iloc[th3]

                # 总坏账率
                Q = (h1 + h2 + df[f"h_{c3}"].iloc[th3]) / 3

                # 此时的最终收入
                temp_value = L * I * P * (1 - Q) - L * P * Q

                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}")

贪心+ QUBO

代码：

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")

————————————————

本文转载自CSDN博主：牵一发而动全身

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

原文链接：https://blog.csdn.net/weixin_41012765/article/details/130249552

一些名词

Ising模型，伊辛模型

QUBO模型，Quadratic Unconstrained Binary Optimization，二次无约束二进制优化模型

1.自旋变量，其域空间取值是Ising空间，也就是域空间是{ − 1 , 1 } 

2.二进制变量，其域空间取值是布尔空间，也就是域空间是{ 0 , 1 } 

HOBO，High Order Binary Optimization，高阶二进制优化问题

组合优化问题，组合优化(Combinatorial Optimization, CO)

二次无约束二进制优化(QUBO)模型

论文参考： Quantum bridge analytics I: a tutorial on formulating and using QUBO models | SpringerLink

可以将很多组合优化问题重新进行表述为QUBO模型，而不同类型的约束关系可以用惩罚函数以非常自然的方式体现在“无约束”QUBO公式中。在QUBO公式中，使用惩罚函数可以产生比较精确的模型表示（这句话不是很懂，翻译过来的）。

QUBO模型用于求解NP难问题，该方法用于设计寻找“最优”解决方案的精确求解器基本不可能，除非是非常小的问题实例。该模型使用现代元启发式方法，在有限的时间内找到接近最优解的方案。

QUBO基本模型

最简单的QUBO模型可以表述为：

其中 x是二进制决策变量的向量，Q 是常数的方阵

这里的常数方阵Q 一般在最优化问题中，在不丧失一般性的情况下，使用对称矩阵或者上/下三角形

对上述进行展开，并改写成标量的形式如下：

由于x _i ∈ { 0 , 1 } 属于二进制变量，存在x _i =，因此对于更具一般性的QUBO模型表示如下:

其中，

QUBO问题可以等价的使用Ising模型来表示，通过变量y_i​=1−2x_i​进行转换，将原本QUBO模型中决策变量域空间{ 0 , 1 }映射到Ising模型决策变量域空间{ − 1 , + 1 }。转变后的Ising模型表述如下：

其中，y_i​∈{−1,1},i∈X

论文中一个小小的案例

最优化问题：最小化二元变量的二次目标函数 y：

其中，变量x_i​∈{0,1}。它有一个线性部分（−5x₁−3x₂−8x₃−6x₄ ）和一个二次部分（4x₁ x₂ +8x₁ x₃ +2x₂ x₃ +10x₃ x₄），由于二进制变量x _j 满足x _j =，所以线性部分可以写成：

基于此，将优化模型的最小化目标函数写成矩阵形式：

由于问题规模较小，穷举(2⁴=16 种可能) 求出该问题的最优解为：y = − 11 , x ₁ = x ₄ = 1 , x ₂ = x ₃ = 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()

这个库有点问题，求的最优解不对，但是可以找到接近最优解。可能是参数设置有问题。

 

使用PyQUBO库：[2103.01708] PyQUBO: Python Library for Mapping Combinatorial Optimization Problems to QUBO Form (arxiv.org)

相关代码：

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模型

一般的QUBO模型要求变量为二进制外不包含任何约束。所以如果要实际解决一些有约束条件的问题，需要在目标函数中引入二次惩罚（quadratic penalties）项。

 这里注意一下P (xy) 表示的是标量P  乘以( x ∗ y )，也就是P ∗ ( x ∗ y )，其他类同。

其中，函数 P计算出来的是正的标量惩罚值，且这个值计算的必须得足够大。因为，在QUBO中，其约束条件是通过优化器来实现的，因此，惩罚项指定的规则是：对于最小化问题的求解，如果是可行解，即满足约束条件，对应的惩罚项等于零；对于不可行解，即不满足约束条件的解，其惩罚项等于一些正的惩罚量。对于P 的取值是否合理可以通过求解模型的解，去计算添加的惩罚函数，如果是0，那就说明，P值效果可以。

惩罚值太大会阻碍求解过程，因为惩罚项会淹没原始目标函数信息，使得很难区分一个解决方案的质量。另一方面，惩罚值过小会危及寻找可行的解决办法。因此，如何确定惩罚值需要进行思考设计。

对于约束项的设计也有相关的技巧，可以查资料看看。

求解

可以将很多组合优化问题重新进行表述为QUBO模型，并且可以使用模拟退火、量子退火和量子近似优化算法(QAOA)等方法对其求解。都是现代元启发式搜索方法。

量子退火

很多专有名词和原理细节可以参考：D-Wave System Documentation或者查资料。比如说纠缠（entanglement）,耦合器（coupler）,本征谱(Eigenspectrum),哈密顿量(Hamiltonian),本征态（eigenstates），初始哈密顿量（Initial Hamiltonian），最终哈密顿量(Final Hamiltonian)，隧穿哈密顿量 (tunneling Hamiltonian)，问题哈密顿量 (problem Hamiltonian)。

量子退火是基于耦合量子位的自然行为来寻找基态（最低能量状态）的过程。量子退火的过程可以使用根据时间变化的哈密顿量来表示：

其中， A ( s ) 和B ( s )是退火路径函数，根据归一化退火时间s = t / t _a ，其中t _a 表示总的退火时间。这样设计可以保证：

那么在设计的时候，初始状态的设计可以根据问题进行简单设计，比如：

其中，表示第 i个量子位的P a u l i − x 算子（泡利− x 算子）

而问题哈密顿量，也就是H _P设计可通过最优化的目标函数给出：

其中，表示第 i个量子位的P a u l i − z 算子（泡利− z算子）。

量子退火器首先初始化量子位的叠加态使得，然后在退火时间内操作量子位耦合，使朝向 H _P 发展，退火结束后，就处于 problem Hamiltonian 的本征态（eigenstate）。根据量子计算的绝热定理，如果退火时间足够长，该时变哈密顿量将始终保持基态，即最优化函数（的解）

模拟退火

模拟退火算法（Simulated Annealing，SA）是一种通用的全局优化算法，用于在搜索空间中找到最优解。它的基本原理是模拟物理学中的退火过程，通过温度参数控制搜索过程中接受次优解的概率，从而避免陷入局部最优解。 模拟退火算法的基本流程如下：

1.初始化一个解 x ₀和一个初始温度 T ₀ 。

2.在该温度下，进行迭代，每次迭代进行以下操作：

a.从当前解 x的邻域中随机选择一个新解 x ′。

b.计算 x ′的目标函数值 f ( x ′ ) 和 x的目标函数值 f ( x )。

c.根据Metropolis准则判断是否接受新的解，即如果 f ( x ′ ) < f ( x )，则以概率1接受 x ′  作为新的解；否则以概率 e − ^{Δ f / T} 接受 x ′ ，其中 Δ f = f ( x ′ ) − f (x) 。

3.更新温度 T，并重复步骤 2 直到温度满足停止条件

其中，温度 T 是一个随着时间逐渐降低的参数，控制着算法接受次优解的概率。

由上述算法可知，其迭代效果受到初始温度、降温速率和终止温度三个参数的影响。

1.初始温度通常设置为一个较高的值，以保证在搜索的早期能够接受一些劣解，从而有更大的概率跳出局部最优解。

2.终止温度通常设置为一个较小的值，当温度降到终止温度时，搜索停止，此时得到的解就是算法得到的最优解。

3.降温速率通常设置为一个小于1的常数，降温速率越慢，则搜索的时间越长，但是搜索范围也就越广，有更大的概率找到全局最优解。

因此需要合理设置初始温度、降温速率和终止温度三个参数以便能够快速达到最优解。

这里介绍一下dwave-samplers库中的SimulatedAnnealingSampler采样器以及对应的参数。

在该库中，初始温度和终止温度是通过转换为β ₁和β ₀ （其中，K是波尔兹常量），通过传入参数 beta_range = [ β 0 , β 1 ] 传入。而降温速率一般通过beta_schedule_type和num_sweeps来设置：对于beta_schedule_type参数，有三种方式：

1. “linear”：线性，在[ β 0 , β 1 ]中线性获取num_sweeps个采样点作为每次更新温度的数值
2. “geometric”，几何，在[ β 0 , β 1 ]中通过np.geomspace(*beta_range, num=num_betas)获取num_sweeps个采样点作为每次更新温度的数值。
3. “custom”：用户自定义

量子近似优化算法(QAOA)

量子近似优化算法可以使用QPanda库，这部分不是很了解，有时间后面补充。

高阶问题（HOBO）

主要解决的问题是：对于QUBO模型中高于二次项的目标函数该如何求解问题

参考论文：[2001.00658] Compressed Quadratization of Higher Order Binary Optimization Problems (arxiv.org)

在不同论文中有不同的叫法：

Higher Order Binary Optimization (HOBO)：高阶二进制优化
termwise quadratizationquadratization:按期限平方，这个是将单项式单次转换为二次。也就是x _i =   , where x _i ∈ 0 , 1 

高阶二元优化问题的压缩二次化：对于目前而言，存在一种使用Rosenberg多项式来减少高阶优化问题次数的方法，该方法在布尔空间中通过引入一个额外的变量来减少一个项的次数，这种方法最终会导致最后的布尔空间矩阵很稀疏，并且问题规模会变大。然后在论文中，提出了一种直接在原有布尔空间中运行的降阶方法，而不通过增加额外变量的方式。

论文中有个没看懂：

感觉这里的x打错了，应该是z，然后想表达的意思应该是，对于z属于布尔空间的取值，可以通过增加变量 y，使得当 y取某个特定值的时候 h ( _{z , y 取特定值} )等价于f ( z )，由于 h函数不确定，只是说能找到该函数满足等价这个条件，那么再做一个限定，也就是当y 取该特定值的时候，正好是h ( _{z , y 取特定值} )取到最小值

那么根据贪心思想，很容易理解，到时候最小化f ( z )等价于最小化 h ( z , y ) 

基于此，就有等价：

这个比较难理解，但是，想到x _i ∈ { 0 , 1 }，并且，这里y是很多个的组合，并且y _i ∈ { 0 , 1 } 就比较容易穷举，在论文Linear and quadratic reformulations of nonlinear optimization problems in binary variables中证明了对于最高次数d ，需要的变量y个数是⌈ log ⁡ d ⌉ − 1。（感觉公式表达有问题）

这种代替方式有两个问题，具体可以看看论文，我们没有使用变量代替降次，所以对于这个HOBO没有仔细研究。

2023年mathorcup的A题

数学定义

问题一

对于问题一，我们定义决策变量 x _{i , j​}用于表示：

那么对于决策变量 x _{i , j}对应的通过率 P _{i , j} 与坏损率 Q _{i , j}有：

于是决策变量 x _{i , j} 对应的最终收入 H _{i , j} 有

考虑整体，也就是最终收入 H有

对于最终收入，我们考虑最小化目标函数 − H ,即

考虑约束条件：由于在整体中只有一个评分卡的阈值能被选取，因此在决策变量x _{i , j} 中只有一个是1，即：

这里的 X空间定义为：X : = { ( i , j ) | i ∈ { 1 , 2 , . . . , 100 } , j ∈ { 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 ，决策出第三张信用评分卡的最佳阈值后，同理决策出其他两张信用卡的最佳阈值。多迭代几轮，这种方法决策出来的阈值在该问题中就是最佳的阈值组合。

具体看下面内容：

QUBO

定义决策变量 x _{i , j , k}用于表示:

x _{i , j , k=}

那么对于决策变量 x _{i , j , k},对应的通过率 P _{i , j , k}与坏损率 Q _{i , j , k}有：

于是决策变量 x _{i , j , k} 对应的最终收入 H _{i , j , k} 有

考虑整体，也就是最终收入 H有

同样考虑最小化目标函数 − H ,此时决策变量的约束条件为：(只有一种方案选取)

这里的 X空间定义为：X : = { ( i , j , k ) |  i , j , k ∈ { 1 , 2 , . . . , 10 } } 

综上所述：

通过添加约束项 λ 将上述二次有约束二进制优化模型改写为二次无约束二进制优化模型（QUBO）:

在代码实现过程中，由于需要计算总通过率T _{1 , i}  ∗  T _{2 , j}   ∗  T _{3 , k}以及总坏账率 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 * df[f"t_{card2}"].iloc[j] * df[f"t_{card3}"].iloc[k]

                # 总坏账率
                Q = (df[f"h_{card1}"].iloc + df[f"h_{card2}"].iloc[j] + df[f"h_{card3}"].iloc[k]) / 3

                # 此时的最终收入
                temp_value = L * I * P * (1 - Q) - L * P * Q

                # 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]")

贪心+ QUBO

在运行过程中，发现一个对第三问很有帮助的规律：可以把上面的决策过程（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]。这里打印出来结果比较少，可能看的不是很明显规律。但是容易发现：

对于三张评分信用卡的决策过程，如果固定第一张和第二张评分信用卡的阈值选取，而考虑第三张评分信用卡的阈值选取，此时选择最优的信用评分卡阈值。然后固定第一张和第三张评分信用卡的阈值选取，而考虑第二张评分信用卡的阈值选取，此时也选择最优的信用评分卡阈值。然后固定第二张和第三张评分信用卡的阈值选取，而考虑第一张评分信用卡的阈值选取，此时也选择最优的信用评分卡阈值。基于此不停的迭代，那么最终获得的阈值选取就是最优的阈值组合。（可以考虑证明一下）。

基于该过程，我们就可以将原本 10 × 10 × 10 个决策变量的求解转换为 迭代次数 × 10 个决策变量的问题。这对于问题三而说是非常重要的一种发现。对于问题三，原本是 1000 × 990 × 980个决策变量问题转换为了 迭代次数 × 980 个决策变量规模的问题。

问题三

最开始我们考虑的是：设计决策变量x _{i , j} 用于表示：

在不考虑约束条件的情况下，整体的通过率P和坏损率Q 满足：

对应的最终收入H有：

可以发现，如果使用这个公式，其目标函数中，决策变量的最高次项系数为1001次。转换为了 HOBO模型，需要降次，感觉很复杂。不考虑。

如果使用这个公式，然后考虑使用多项式近似(Polynomial approximation)的思想，将指数项使用多项式展开，比如泰勒级数展开式：,保留一阶项，即：。也就可以将指数项转换为一个多项式函数的线性项：

那么基于这个理论就可以将原本的目标函数转换为最高二次项的目标函数。不过这个也有点问题，就是由于决策变量x _{i , j} ∈ { 0 , 1 } 是离散变量，那么使用泰勒级展开式未必等价， 也不做考虑。

以上所有考虑都存在一个问题，就是怎样解决整体的通过率P表述。如果能够将P 转换为一次项的表达，那么对于整个求解问题就很简单了。那么就重新定义决策变量x _{i , j , k , l , m , n}用于决策是否选择：第 i张信用评分卡选第 l 个阈值，第 j张信用评分卡选第 m 个阈值，第 k 张信用评分卡选第 n个阈值。即：

那么对于决策变量 x _{i , j , k , l , m , n} ，对应的通过率 P _{i , j , k , l , m , n} 与坏损率 Q _{i , j , k , l , m , n} 有：

那么对于决策变量 x _{i , j , k , l , m , n} 对应的通过率 P _{i , j , k , l , m , n} 与坏损率 Q _{i , j , k , l , m , n} 有：

那么代入最终收入公式，有：

=

这种方法如果不考虑约束条件，那么量子个数有 100 × 100 × 100 × 10 × 10 × 10 ，内存有点不太够，运行不出来这个。

考虑暴力

考虑使用循环暴力遍历其找到最佳的组合。这个需要很长的时间来进行，那么考虑多线程+分治的思想对其暴力，代码如下：

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):
                            # 总通过率
                            P = df[f"t_{i}"].iloc[m] * df[f"t_{j}"].iloc[n] * df[f"t_{k}"].iloc[l]

                            # 总坏账率
                            Q = (df[f"h_{i}"].iloc[m] + df[f"h_{j}"].iloc[n] + df[f"h_{k}"].iloc[l]) / 3

                            # 此时的最终收入
                            temp_value = L * I * P * (1 - Q) - L * P * Q

                            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):
                # 总通过率
                P = t1 * t2 * df[f"t_{c3}"].iloc[th3]

                # 总坏账率
                Q = (h1 + h2 + df[f"h_{c3}"].iloc[th3]) / 3

                # 此时的最终收入
                temp_value = L * I * P * (1 - Q) - L * P * Q

                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}")

贪心+ QUBO

代码：

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")

————————————————

本文转载自CSDN博主：牵一发而动全身

版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。

原文链接：https://blog.csdn.net/weixin_41012765/article/details/130249552