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