数学建模|通过模拟退火算法求解供应与选址问题:问题一(python代码实现)

2024-01-07 17:41:30

今天继续来学习模拟退火算法在数学建模中的应用,如果对模拟退火算法的基础知识还不了解的,可以看我之前的博客。

通过模拟退火算法求解一元五次方程最值(python代码实现)-CSDN博客

这次要解决的供应与选址问题依然来自数学建模老哥的视频:

13 非线性规划算法在数学建模中的应用与编程实现_哔哩哔哩_bilibili

问题如下:

如果对这个问题还不是很了解,可以先去看视频,我在这里就不过多解释。我在这里主要解决用编程求解这个问题。

首先看到第一问。先把这个问题转化为一个规划问题,求一个最小值。那么,视频里已经帮我们转化好了,如下:

别看他写的那么复杂,其实目标函数就是距离乘供货量,这里画了个图,可以感受一下(画的不好~)。现在AB的坐标和六个工地的坐标都是固定的,所以距离也是固定的。现在的问题就变成了,A和B要向这六个工地供货多少的问题了。

通过上面的分析,我们知道了,A和B的供货量,就是我们要求的自变量。我们可以先用一个numpy数组存储一下我们的解,也就是自变量。

这里生成了一个两行六列的全0数组。第一行表示A向六个工地的供货量,第二行表示B向六个工地的供货量。

我们再用numpy数组计算一下A和B到六个工地的距离。

这里还是用两个numpy数组存储了工地和料场的信息,第一行和第二行是工地和料场的坐标,gondi的第三行,是各个工地的水泥需求量。

这里numpy的函数就能计算出距离了,结果如下:

现在距离有了,自变量有了,我们的目标函数就可以写出来了。

target_fun = lambda x,juli: np.sum(juli * x) #注意这里是numpy的数组乘法,不是矩阵乘法

目标函数有了,剩下的就是约束条件。这里自变量受到两个约束:1.A和B的供货量不能超过自身的储备量,也就是不能超过20。2.A和B的供货量要等于一个工地自身的需求量。

反应到自变量数组上就是,数组的每一行的所有元素相加要小于等于20,数组的每一列的元素相加要刚好等于对应工地的需求量。

根据我们的分析,写出我们的约束条件,这里用numpy计算和判断就很方便了。

st_fun1 = lambda x: np.sum(x,axis=1)[0]<=20 and np.sum(x,axis=1)[1]<=20
st_fun2 = lambda x ,gondi: np.all(np.sum(x,axis=0) == gondi[2], axis=0)
judgy_st = lambda x,gondi:st_fun1(x) and st_fun2(x,gondi)

接着,我们来设计自变量的获取方式和更新方式。这里有很多的坑,我们来一一分析。

首先,我想到的就是,先随机取数组的一个地址,然后给那个地址赋予一个随机浮点数。类似于这样(补充一点,还要注意自变量的上下限,下限是0,上限是工地对应的需求量):

def get_x(x,lb,t_move):
    i = random.randint(0,1)
    j = random.randint(0,5)
    ub = gondi[2][j]
    lb = x[i][j]-t_move*(x[i][j]-lb)
    ub = x[i][j]+t_move*(ub-x[i][j])
    # print(t_move)
    x[i][j] = random.uniform(lb,ub)
    return x

但这样是不行的,这样很难得到我们想要的符合约束条件的解。为了优化这个问题,我进一步想到的了,每一列的元素之间是有联系的:列元素相加就等于对应工地的需求量。也就是说,为某一列的一个元素赋值,那么这一列的另外一个元素的值也就出来了,等于 需求量 - 随机赋的值。于是我把代码修改成这样:

def get_x(x,lb,t_move):
    i = random.randint(0,1)
    j = random.randint(0,5)
    ub = gondi[2][j]
    lb = x[i][j]-t_move*(x[i][j]-lb)
    ub = x[i][j]+t_move*(ub-x[i][j])
    # print(t_move)
    x[i][j] = random.uniform(lb,ub)
    x[(i+1) % 2][j] = ub-x[i][j]
    return x

这样就能很快得到我们想要的解了。但是,还有问题。就是我们得到的解的元素都是浮点数,很难得到整数解。通过视频开天眼,我们可以知道最后的解的元素都是整数。通过浮点数求出的解并不是最优解。我继续修改代码,让它能达到一个随机整数:

def get_x(x,lb,t_move):
    i = random.randint(0,1)
    j = random.randint(0,5)
    ub = gondi[2][j]
    lb = x[i][j]-t_move*(x[i][j]-lb)
    ub = x[i][j]+t_move*(ub-x[i][j])
    # print(t_move)
    if random.random() > random.random():
        x[i][j] = random.randint(int(lb), int(ub))
        x[(i + 1) % 2][j] = ub - x[i][j]
    else:
        x[i][j] = random.uniform(lb,ub)
        x[(i+1) % 2][j] = ub-x[i][j]
    return x

这就是最终的get_x函数啦,通过这样的设置,它就有一半的概率得到一个整数解,有一半的概率得到一个浮点数解。

完成了定义函数和设计自变量的获取方式,我们就可以直接带入模拟退火算法进行求解了。至于怎么带入,可以看我之前的博客。

数学建模|通过模拟退火算法求解非线性规划问题(python代码实现)-CSDN博客

这里直接给出完整代码:

#%%
import random
import numpy as np
import math
t = 200000000
T = 200000000
dt = 0.999999993
eps = 1e-14
n = 0

gondi = np.array([
    [1.25,8.75,0.5,5.75,3,7.25],
    [1.25,0.75,4.75,5,6.5,7.25],
    [3,5,4,7,6,11]
])

liaochang = np.array([
    [5,2],
    [1,7],
])

#%%

juli = np.array([
    np.sqrt(np.sum((gondi[0:2,:]-liaochang[0:2,0:1])**2,axis=0)),
    np.sqrt(np.sum((gondi[0:2,:]-liaochang[0:2,1:2])**2,axis=0))
])

#%%

x = np.zeros([2,6])

target_fun = lambda x,juli: np.sum(juli * x)
st_fun1 = lambda x: np.sum(x,axis=1)[0]<=20 and np.sum(x,axis=1)[1]<=20
st_fun2 = lambda x ,gondi: np.all(np.sum(x,axis=0) == gondi[2], axis=0)
judgy_st = lambda x,gondi:st_fun1(x) and st_fun2(x,gondi)


#%%
def get_x(x,lb,t_move):
    i = random.randint(0,1)
    j = random.randint(0,5)
    ub = gondi[2][j]
    lb = x[i][j]-t_move*(x[i][j]-lb)
    ub = x[i][j]+t_move*(ub-x[i][j])
    # print(t_move)
    if random.random() > random.random():
        x[i][j] = random.randint(int(lb), int(ub))
        x[(i + 1) % 2][j] = ub - x[i][j]
    else:
        x[i][j] = random.uniform(lb,ub)
        x[(i+1) % 2][j] = ub-x[i][j]
    return x

#%%
move = 11
jud_n = 0
best_x = get_x(x,0,1)
# breakpoint()
while judgy_st(best_x,gondi) != True:
    best_x = get_x(best_x,0,1)
    jud_n += 1
best_y = target_fun(best_x,juli)
print(f"f = {best_y},x = {best_x},第{jud_n}次循环")
y = best_y

while n!= 200000 :
    # print(t/T)
    dx = get_x(x,0,1)
    judgy_n = 0
    while judgy_st(dx,gondi) != True:
        dx = get_x(dx,0,1)
        judgy_n += 1
        if judgy_n == 1000000:
            break
    if judgy_n == 1000000:
        continue
    dy = target_fun(dx,juli)

    if dy < y:
        y = dy
        x = dx
        if dy < best_y:
            best_y = dy
            best_x = dx
            print(f"f = {best_y},\n x = {best_x},\n第{n}次循环")

    elif  np.exp((y - dy) / t)  > random.random():
        y = dy
        x = dx

    # print(f"f = {y},x = {x},第{n}次循环")

    t = t * dt
    n += 1

print(n)

最后我们来看一下结果:

这是视频的计算结果:

可以看到,我们的结果与视频的结果还是比较相近的,甚至我们计算的还要更小一点。

想看问题二的编程思路的,可以看这里。

数学建模|通过模拟退火算法求解供货与选址问题:问题二(python代码实现)-CSDN博客

文章来源:https://blog.csdn.net/2302_80902795/article/details/135407921
本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。