Q1: 单信用评分卡求最大收入

符号定义

C：总贷款资金

R：利息收入率

P：惩罚系数（约束条件的）

:二进制变量，表示是否选择第 i 张信用评分卡

:表示是否选择第 i 张信用评分卡的第 j 个阈值

:表示第 i 张信用评分卡的第 j 个阈值对应的通过率

:表示第 i 张信用评分卡的第 j 个阈值对应的坏账率

建模

1.定义决策变量：

用一个二进制变量表示是否选择第 i 张信用评分卡，如果选择则=1，否则=0。我们还可以用一个二进制变量表示是否选择第 i 张信用评分卡的第 j 个阈值，如果选择则=1，否则=0。这样我们就有 100 + 1000 = 1100 个决策变量。

2.定义目标函数：

我们的目标是最大化银行的最终收入，根据题目给出的公式，我们可以写出目标函数为：

maxM=，其中和分别表示第 i 张信用评分卡的第 j 个阈值对应的通过率和坏账率，贷款资金和利息收入率是已知的常数。

3.定义约束条件：

需要满足以下约束条件：

(1) 只能选择一张信用评分卡，即

(2) 每张信用评分卡只能选择一个阈值，即对于任意 i，有

(3) 所有的决策变量都是二进制的，即对于任意 i 和 j，有∈ 0 , 1和∈ 0 , 1

忽略常数项 C:

还可以进行一些简化:

因为仅选择一张卡，所以目标函数中的是可以省去的，也就是说只使用来表示是否选择第i张卡的第j个阈值，此时决策变量从1100变为了1000。如果=1，那么表示选择了第i张卡的第j个阈值，现在我们可以将约束条件整合成就行。

按以下规则拼接为，以卡号的顺序排列:

则有

将i,j也表示为k的形式:

QUBO形式

为了将目标函数化为 QUBO 形式，需要将所有的一次项和常数项都转化为二次项。可以利用以下的等式来实现这一目的：

1.=，因为是二值变量，所以它的平方等于它本身。

2.

将目标函数化为 ：

最后将其转化为最小化，忽略常数项P:

此时可以用QUBO矩阵Q表示，有:

其中 Q的元素为：

结果

最终的结果是选择第 49 张卡的第 1 个阈值，对应的收入为 61172.0

代码

这一部分的求解有很多种方式，比如：你可以使用 D-wave 提供的 API 在量子计算机上真正的进行求解，也可以使用OpenJij库的函数 SQASampler() 模拟量子退火求解。

1. 导入需要用到的库

# 数据处理
import numpy as np
import pandas as pd

# 导入 PyQUBO 库
from pyqubo import Array, Binary, Constraint

# 模型求解
from dwave.system import LeapHybridSampler
from openjij import SQASampler

2. 数据处理

# 读取信用卡数据
data = pd.read_csv("附件1：data_100.csv", header=0)
t = data.iloc[:, ::2].values.T   # 通过率
h = data.iloc[:, 1::2].values.T  # 坏账率
print(t.shape, h.shape)          # 行对应阈值，列对应卡

3. 定义函数处理后面的结果

import re

def ret_result(solution):
    '''
    返回索引和结果
    '''
    i,j = 0,0
    pattern_idx = re.compile(r"[0-9]+")
    for key,value in solution.items():
        if value == 1:
            print(key, value)
            idx = pattern_idx.findall(key)
            if len(idx)==1:
                i = int(idx[0])
                print(i)

            else:
                i, j = int(idx[0]), int(idx[1])
                print(i, j)

    result = C * R * t[i, j] * (1 - h[i, j]) - C * t[i, j] * h[i, j]
    return i, j, result

4. 模型参数设置和求解

# 定义常数和变量
C = 1000000 # 总贷款资金
R = 0.08    # 利率
n = 100     # 信用评分卡数量
P = C       # 惩罚系数

# 定义目标函数
H = 0
for i in range(100):
    for j in range(10):
        H -= x * y[i, j] * (C * R * t[i, j] * (1 - h[i, j]) - C * t[i, j] * h[i, j])

# 定义约束条件
C1 = P * (sum(x) - 1)**2
C2 = P * sum((sum(y) - x)**2 for i in range(100))

# 将目标函数和约束条件相加
model = (H + C1 + C2).compile()

# 利用to_qubo()得到矩阵Q
Q, _ = model.to_qubo()

# 求解 QUBO 模型
sampler = LeapHybridSampler(num_reads=3)        # 使用 LeapHybridSampler
response = sampler.sample_qubo(Q) # 采样 QUBO
solution = response.first.sample     # 获取最优解
i, j, result = ret_result(solution)
print(f"第{i+1}张卡选择第{j+1}个阈值的收入是: {result}")

5. 转化为单一变量后的代码

z = [Binary(f"z_{k}") for k in range(n*10+100)]

# 定义目标函数
H = 0
for i in range(100):
    for j in range(10):
        H -= z * z[10 * i + j + 100] * (C * R * t[i, j] * (1 - h[i, j]) - C * t[i, j] * h[i, j])

# 定义约束条件
C1 = P * (sum(z[:100]) - 1)**2
C2 = P * sum((sum(z[i * 10 + 100: i*10+110]) - z)**2 for i in range(100))

# 将目标函数和约束条件相加
model = (H + C1 + C2).compile()

# 利用to_qubo()得到矩阵Q
Q, _ = model.to_qubo()

# 求解 QUBO 模型
sampler = LeapHybridSampler()        # 使用 LeapHybridSampler
response = sampler.sample_qubo(Q)    # 采样 QUBO
solution = response.first.sample     # 获取最优解

# 打印对应的索引
pattern_idx = re.compile(r"[0-9]+")
for key,value in solution.items():
    if value == 1:
        print(key, value)
        idx = int(pattern_idx.findall(key)[0])
        print(idx)

 

Q2：前三张信用评分卡求解最大收入

建模

因为已经给定了三张信用评分卡，只需要选择每张信用评分卡的阈值。

可以用三个二值向量，，来表示三张信用评分卡的阈值，每个变量有 10 种取值，分别对应 10 个阈值。

目标函数是最大化银行收入减去坏账损失，约束条件是每个变量只能取一个值。

1.定义决策变量：

为一个二值变量，表示是否选择第 i 张信用评分卡的第 j 个阈值，其中 i = 1, 2, 3，j = 1, 2, …, 10。

2.定义目标函数：

我们的目标是最大化银行的最终收入，根据题目给出的公式，我们可以写出目标函数为：

因为该题规定了选择前三张卡，所以我们可以简化总通过率和总坏账率的表达形式为:

忽略常数项做简化，有:

3.定义约束条件：

每张信用评分卡只能选择一个阈值，即：其中 i = 1, 2, 3。

将约束条件加入目标函数中，并引入一个惩罚项系数 P，得到：

QUBO的形式–三次转二次

这个目标函数是三次的，我们需要将其转为二次，这里设*=，这一操作对应的约束条件为：

这些约束条件可以保证当且仅当=1且=1时，=1。

根据参考文献（A Tutorial on Formulating and Using QUBO Models）中的

可以得到对应的惩罚项：

汇总即有：

将惩罚项加入目标函数：

打断一下，下面的拼接只是为了写出 QUBO 矩阵，其实能弄明白第一题的怎么写就行，不需要太花费时间，这里只是一个简单的二次型矩阵表达，但因为拼接后看起来会很复杂。

现在，我们需要将，，和合并成一个向量，拼接规则如下:

此时，对应的约束条件转化为：

目标函数：

因为

可以将目标函数转换为以下形式：

然后，可以将上面的表达式转换为 QUBO 矩阵的形式。具体来说，可以定义一个 130 × 130的矩阵 Q ，其中 Q 的元素为：

结果

最终的结果是分别选择前三张卡的第 8，1，2 个阈值，最终的收入为27914.8。

代码

1. 导入需要的库以及数据处理

# 导入需要的库
import re

import numpy as np
import pandas as pd

from openjij import SQASampler
# from dwave.system import DWaveSampler, LeapHybridSampler, EmbeddingComposite
from pyqubo import Array, Binary


# 读取信用卡数据
data = pd.read_csv("附件1：data_100.csv", header=0)
t = data.iloc[:, ::2].values.T   # 通过率
h = data.iloc[:, 1::2].values.T  # 坏账率
print(t.shape, h.shape)          # 行对应卡，列对应阈值

2. 定义函数处理结果

def get_result(solution):
    '''
    返回结果并打印索引及结果
    '''
    # x, y 用于存储下标
    x, y = [], []
    i = 0
    pattern_idx = re.compile(r"[0-9]+")
    for key, value in solution.items():
        if len(key) < 10:
            if value == 1 and i < 3:
                print(key, value)
                x_idx, y_idx = pattern_idx.findall(key)
                print(f"第{int(x_idx)+1}张卡选择第{int(y_idx)+1}个阈值")
                x.append(int(x_idx))
                y.append(int(y_idx))
                i += 1
    # T, H 分别表示总通过率和总坏账率
    T = t[x[0], y[0]] * t[x[1], y[1]] * t[x[2], y[2]]
    H = 1 / 3 * (h[x[0], y[0]] + h[x[1], y[1]] + h[x[2], y[2]])
    result = C * T * (R * (1 - H) - H)
    print(f"result\n")
    return result

3. 求解

# 定义常数和变量
C = 1000000  # 总贷款资金
R = 0.08     # 利率
n = 100      # 信用评分卡数量
P = C        # 惩罚系数

# 定义目标函数
y = Array.create('y', shape=(3, 10), vartype='BINARY')
z = Array.create('z', shape=(10, 10), vartype='BINARY')

Q = 0
for i in range(10):
    for j in range(10):
        for k in range(10):
            T_ijk = t[0, i] * t[1, j] * t[2, k]
            H_ijk = (h[0, i] + h[1, j] + h[2, k]) / 3
            Q += y[0, i] * z[j, k] * T_ijk * (C * R * (1 - H_ijk) - C * H_ijk)

c1 = P * sum((sum(y) - 1)**2 for i in range(3))
c2 = 0
for j in range(10):
    for k in range(10):
        c2 += P * (3 * z[j, k] - 2 * z[j, k] * y[1, j] -
                   2 * z[j, k] * y[2, k] + y[1, j] * y[2, k])

model = -Q + c1 + c2
model = model.compile()

qubo, offset = model.to_qubo()

best_result = 0
results = []
print("因为输出不稳定，所以写了个循环代码进行结果的验证，你可以修改SQASampler来得到更稳定的结果（每次运行大概需要一分半）\n")
for k in range(20):
    # sampler = LeapHybridSampler()  # 使用 LeapHybridSampler，如果没有api的话就用下面那行
    sampler = SQASampler(
        num_sweeps=3000,
        num_reads=5000)  # num_sweeps每个温度的扫描次数，num_reads指重复模拟退火算法的次数
    response = sampler.sample_qubo(qubo, )  # 采样 QUBO
    solution = response.first.sample
    
    print(f"第{k+1}次运行结果：")
    result = get_result(solution)
    results.append(result)

    if result > best_result:
        best_result = result
        best_solution = solution

print("最优解:")
get_result(best_solution)

# 单次运行结果（注释上面循环部分，再取消注释下面）
# sampler = SQASampler(num_sweeps=3000, num_reads=5000)  # num_sweeps每个温度的扫描次数，num_reads指重复模拟退火算法的次数
# response = sampler.sample_qubo(qubo, )  # 采样 QUBO
# solution = response.first.sample

# get_result(solution)

需要注意的是，代码中模拟量子退火的解是不稳定的，在我最近一次的试运行中，10 次中输出最优解的次数为 7 次。
————————————————

本文转载自CSDN博主：Hoper. J

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

原文链接：https://blog.csdn.net/weixin_42426841/article/details/130209098

Q1: 单信用评分卡求最大收入

符号定义

C：总贷款资金

R：利息收入率

P：惩罚系数（约束条件的）

:二进制变量，表示是否选择第 i 张信用评分卡

:表示是否选择第 i 张信用评分卡的第 j 个阈值

:表示第 i 张信用评分卡的第 j 个阈值对应的通过率

:表示第 i 张信用评分卡的第 j 个阈值对应的坏账率

建模

1.定义决策变量：

用一个二进制变量表示是否选择第 i 张信用评分卡，如果选择则=1，否则=0。我们还可以用一个二进制变量表示是否选择第 i 张信用评分卡的第 j 个阈值，如果选择则=1，否则=0。这样我们就有 100 + 1000 = 1100 个决策变量。

2.定义目标函数：

我们的目标是最大化银行的最终收入，根据题目给出的公式，我们可以写出目标函数为：

maxM=，其中和分别表示第 i 张信用评分卡的第 j 个阈值对应的通过率和坏账率，贷款资金和利息收入率是已知的常数。

3.定义约束条件：

需要满足以下约束条件：

(1) 只能选择一张信用评分卡，即

(2) 每张信用评分卡只能选择一个阈值，即对于任意 i，有

(3) 所有的决策变量都是二进制的，即对于任意 i 和 j，有∈ 0 , 1和∈ 0 , 1

忽略常数项 C:

还可以进行一些简化:

因为仅选择一张卡，所以目标函数中的是可以省去的，也就是说只使用来表示是否选择第i张卡的第j个阈值，此时决策变量从1100变为了1000。如果=1，那么表示选择了第i张卡的第j个阈值，现在我们可以将约束条件整合成就行。

按以下规则拼接为，以卡号的顺序排列:

则有

将i,j也表示为k的形式:

QUBO形式

为了将目标函数化为 QUBO 形式，需要将所有的一次项和常数项都转化为二次项。可以利用以下的等式来实现这一目的：

1.=，因为是二值变量，所以它的平方等于它本身。

2.

将目标函数化为 ：

最后将其转化为最小化，忽略常数项P:

此时可以用QUBO矩阵Q表示，有:

其中 Q的元素为：

结果

最终的结果是选择第 49 张卡的第 1 个阈值，对应的收入为 61172.0

代码

这一部分的求解有很多种方式，比如：你可以使用 D-wave 提供的 API 在量子计算机上真正的进行求解，也可以使用OpenJij库的函数 SQASampler() 模拟量子退火求解。

1. 导入需要用到的库

# 数据处理
import numpy as np
import pandas as pd

# 导入 PyQUBO 库
from pyqubo import Array, Binary, Constraint

# 模型求解
from dwave.system import LeapHybridSampler
from openjij import SQASampler

2. 数据处理

# 读取信用卡数据
data = pd.read_csv("附件1：data_100.csv", header=0)
t = data.iloc[:, ::2].values.T   # 通过率
h = data.iloc[:, 1::2].values.T  # 坏账率
print(t.shape, h.shape)          # 行对应阈值，列对应卡

3. 定义函数处理后面的结果

import re

def ret_result(solution):
    '''
    返回索引和结果
    '''
    i,j = 0,0
    pattern_idx = re.compile(r"[0-9]+")
    for key,value in solution.items():
        if value == 1:
            print(key, value)
            idx = pattern_idx.findall(key)
            if len(idx)==1:
                i = int(idx[0])
                print(i)

            else:
                i, j = int(idx[0]), int(idx[1])
                print(i, j)

    result = C * R * t[i, j] * (1 - h[i, j]) - C * t[i, j] * h[i, j]
    return i, j, result

4. 模型参数设置和求解

# 定义常数和变量
C = 1000000 # 总贷款资金
R = 0.08    # 利率
n = 100     # 信用评分卡数量
P = C       # 惩罚系数

# 定义目标函数
H = 0
for i in range(100):
    for j in range(10):
        H -= x * y[i, j] * (C * R * t[i, j] * (1 - h[i, j]) - C * t[i, j] * h[i, j])

# 定义约束条件
C1 = P * (sum(x) - 1)**2
C2 = P * sum((sum(y) - x)**2 for i in range(100))

# 将目标函数和约束条件相加
model = (H + C1 + C2).compile()

# 利用to_qubo()得到矩阵Q
Q, _ = model.to_qubo()

# 求解 QUBO 模型
sampler = LeapHybridSampler(num_reads=3)        # 使用 LeapHybridSampler
response = sampler.sample_qubo(Q) # 采样 QUBO
solution = response.first.sample     # 获取最优解
i, j, result = ret_result(solution)
print(f"第{i+1}张卡选择第{j+1}个阈值的收入是: {result}")

5. 转化为单一变量后的代码

z = [Binary(f"z_{k}") for k in range(n*10+100)]

# 定义目标函数
H = 0
for i in range(100):
    for j in range(10):
        H -= z * z[10 * i + j + 100] * (C * R * t[i, j] * (1 - h[i, j]) - C * t[i, j] * h[i, j])

# 定义约束条件
C1 = P * (sum(z[:100]) - 1)**2
C2 = P * sum((sum(z[i * 10 + 100: i*10+110]) - z)**2 for i in range(100))

# 将目标函数和约束条件相加
model = (H + C1 + C2).compile()

# 利用to_qubo()得到矩阵Q
Q, _ = model.to_qubo()

# 求解 QUBO 模型
sampler = LeapHybridSampler()        # 使用 LeapHybridSampler
response = sampler.sample_qubo(Q)    # 采样 QUBO
solution = response.first.sample     # 获取最优解

# 打印对应的索引
pattern_idx = re.compile(r"[0-9]+")
for key,value in solution.items():
    if value == 1:
        print(key, value)
        idx = int(pattern_idx.findall(key)[0])
        print(idx)

 

Q2：前三张信用评分卡求解最大收入

建模

因为已经给定了三张信用评分卡，只需要选择每张信用评分卡的阈值。

可以用三个二值向量，，来表示三张信用评分卡的阈值，每个变量有 10 种取值，分别对应 10 个阈值。

目标函数是最大化银行收入减去坏账损失，约束条件是每个变量只能取一个值。

1.定义决策变量：

为一个二值变量，表示是否选择第 i 张信用评分卡的第 j 个阈值，其中 i = 1, 2, 3，j = 1, 2, …, 10。

2.定义目标函数：

我们的目标是最大化银行的最终收入，根据题目给出的公式，我们可以写出目标函数为：

因为该题规定了选择前三张卡，所以我们可以简化总通过率和总坏账率的表达形式为:

忽略常数项做简化，有:

3.定义约束条件：

每张信用评分卡只能选择一个阈值，即：其中 i = 1, 2, 3。

将约束条件加入目标函数中，并引入一个惩罚项系数 P，得到：

QUBO的形式–三次转二次

这个目标函数是三次的，我们需要将其转为二次，这里设*=，这一操作对应的约束条件为：

这些约束条件可以保证当且仅当=1且=1时，=1。

根据参考文献（A Tutorial on Formulating and Using QUBO Models）中的

可以得到对应的惩罚项：

汇总即有：

将惩罚项加入目标函数：

打断一下，下面的拼接只是为了写出 QUBO 矩阵，其实能弄明白第一题的怎么写就行，不需要太花费时间，这里只是一个简单的二次型矩阵表达，但因为拼接后看起来会很复杂。

现在，我们需要将，，和合并成一个向量，拼接规则如下:

此时，对应的约束条件转化为：

目标函数：

因为

可以将目标函数转换为以下形式：

然后，可以将上面的表达式转换为 QUBO 矩阵的形式。具体来说，可以定义一个 130 × 130的矩阵 Q ，其中 Q 的元素为：

结果

最终的结果是分别选择前三张卡的第 8，1，2 个阈值，最终的收入为27914.8。

代码

1. 导入需要的库以及数据处理

# 导入需要的库
import re

import numpy as np
import pandas as pd

from openjij import SQASampler
# from dwave.system import DWaveSampler, LeapHybridSampler, EmbeddingComposite
from pyqubo import Array, Binary


# 读取信用卡数据
data = pd.read_csv("附件1：data_100.csv", header=0)
t = data.iloc[:, ::2].values.T   # 通过率
h = data.iloc[:, 1::2].values.T  # 坏账率
print(t.shape, h.shape)          # 行对应卡，列对应阈值

2. 定义函数处理结果

def get_result(solution):
    '''
    返回结果并打印索引及结果
    '''
    # x, y 用于存储下标
    x, y = [], []
    i = 0
    pattern_idx = re.compile(r"[0-9]+")
    for key, value in solution.items():
        if len(key) < 10:
            if value == 1 and i < 3:
                print(key, value)
                x_idx, y_idx = pattern_idx.findall(key)
                print(f"第{int(x_idx)+1}张卡选择第{int(y_idx)+1}个阈值")
                x.append(int(x_idx))
                y.append(int(y_idx))
                i += 1
    # T, H 分别表示总通过率和总坏账率
    T = t[x[0], y[0]] * t[x[1], y[1]] * t[x[2], y[2]]
    H = 1 / 3 * (h[x[0], y[0]] + h[x[1], y[1]] + h[x[2], y[2]])
    result = C * T * (R * (1 - H) - H)
    print(f"result\n")
    return result

3. 求解

# 定义常数和变量
C = 1000000  # 总贷款资金
R = 0.08     # 利率
n = 100      # 信用评分卡数量
P = C        # 惩罚系数

# 定义目标函数
y = Array.create('y', shape=(3, 10), vartype='BINARY')
z = Array.create('z', shape=(10, 10), vartype='BINARY')

Q = 0
for i in range(10):
    for j in range(10):
        for k in range(10):
            T_ijk = t[0, i] * t[1, j] * t[2, k]
            H_ijk = (h[0, i] + h[1, j] + h[2, k]) / 3
            Q += y[0, i] * z[j, k] * T_ijk * (C * R * (1 - H_ijk) - C * H_ijk)

c1 = P * sum((sum(y) - 1)**2 for i in range(3))
c2 = 0
for j in range(10):
    for k in range(10):
        c2 += P * (3 * z[j, k] - 2 * z[j, k] * y[1, j] -
                   2 * z[j, k] * y[2, k] + y[1, j] * y[2, k])

model = -Q + c1 + c2
model = model.compile()

qubo, offset = model.to_qubo()

best_result = 0
results = []
print("因为输出不稳定，所以写了个循环代码进行结果的验证，你可以修改SQASampler来得到更稳定的结果（每次运行大概需要一分半）\n")
for k in range(20):
    # sampler = LeapHybridSampler()  # 使用 LeapHybridSampler，如果没有api的话就用下面那行
    sampler = SQASampler(
        num_sweeps=3000,
        num_reads=5000)  # num_sweeps每个温度的扫描次数，num_reads指重复模拟退火算法的次数
    response = sampler.sample_qubo(qubo, )  # 采样 QUBO
    solution = response.first.sample
    
    print(f"第{k+1}次运行结果：")
    result = get_result(solution)
    results.append(result)

    if result > best_result:
        best_result = result
        best_solution = solution

print("最优解:")
get_result(best_solution)

# 单次运行结果（注释上面循环部分，再取消注释下面）
# sampler = SQASampler(num_sweeps=3000, num_reads=5000)  # num_sweeps每个温度的扫描次数，num_reads指重复模拟退火算法的次数
# response = sampler.sample_qubo(qubo, )  # 采样 QUBO
# solution = response.first.sample

# get_result(solution)

需要注意的是，代码中模拟量子退火的解是不稳定的，在我最近一次的试运行中，10 次中输出最优解的次数为 7 次。
————————————————

本文转载自CSDN博主：Hoper. J

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

原文链接：https://blog.csdn.net/weixin_42426841/article/details/130209098