数学建模|用python求解非线性规划问题（模拟退火算法实现）

对模拟退火基础知识不了解的，可以先去看这个视频，我觉得讲的还是比较通俗易懂的

大学生速通模拟退火算法_哔哩哔哩_bilibili

这次要解决的问题来自数学建模老哥的视频

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

题目如下：

通过模拟退火算法用python求解这个问题，有以下几个步骤：

1. 搭建模拟退火算法的框架
2. 将数学公式转写为python函数
3. 设计自变量的获取方式和更新方式
4. 带入算法进行计算

搭建模拟退火算法的框架

import random
import numpy as np
import math
t = 2000 #t是当前温度
T = 2000 #T是初始温度
dt = 0.993 #dt是温度的变化率
eps = 1e-14
n = 0 #n是外循环的运行次数
x = [0,0,0] #初始化解，用列表存储


while n < 200 : #这里的数字可以设置为希望外循环运行的次数
    #中间是退火代码的实现方式，先留空

    t = t * dt
    n += 1

将数学公式转写为python函数

定义目标函数：

def target_fun(x:list):
    return x[0]**2 + x[1] ** 2 + x[2]**2 + 8

定义约束条件：

def st_fun1(x:list):
    if x[0]**2 - x[1] + x[2]**2 >= 0:
        return True
    else:
        return False

这里可以将函数定义成lambda表达式，可以让我们的代码代码更简洁，这样写的效果和上面写的效果和用法是一样的：

target_fun = lambda x: x[0]**2 + x[1] ** 2 + x[2]**2 + 8
st_fun1 = lambda x: x[0]**2 - x[1] + x[2]**2 >= 0
st_fun2 = lambda x:  x[0] + x[1]**2 +x[2]**2 <= 20
st_fun3 = lambda x:  math.isclose(-x[0] - x[1]**2 + 2 , 0,rel_tol=0.1, abs_tol=0.001)
st_fun4 = lambda x:  math.isclose(x[1] + 2*x[2]**2 , 3,rel_tol=1e-1, abs_tol=0.01)
st_fun5 = lambda x: x[0] >= 0 and x[1] >= 0 and x[2] >= 0
judgy_st = lambda x: st_fun1(x) and st_fun2(x) and st_fun3(x) and st_fun4(x) and st_fun5(x)

注意第三个约束条件和第四个约束条件是一个等式，但在计算过程中是很难严格相等的，因此要用math.isclose()在误差允许的范围内近似相等，就相当于约等于。math.isclose()函数的功能如下：

math.isclose(a, b, rel_tol=1e-9, abs_tol=0.0)

其中，a和b是要比较的两个数值，rel_tol是相对容差（默认值为1e-9），abs_tol是绝对容差（默认值为0）。如果两个数值的差小于等于这两个容差的任意一个，则认为它们近似相等。

当 a-b<rel_tol 或者 |a-b|<abs_tol 时就可认为他们相等

最后还要添加一个judge_st()函数，判断我们的解是否满足所有的约束条件。

设计自变量的获取方式和更新方式

代码如下：

move = 5 #x的定义域
t_move = t/T #当前温度除以初始温度的比率，随着温度降低而减小

def get_x(x:list,lb,ub,t_move): #随机获取一个浮点数作为解
    i = random.randint(0,2)
    lb = x[i]-t_move*(x[i]-lb)
    ub = x[i]+t_move*(ub-x[1])
    # print(t_move)
    x[i] = random.uniform(lb,ub)
    return x

jud_n = 0
best_x = get_x(x,0,move,t/T) #这里先初始化一个最优解
# breakpoint()
while judgy_st(best_x) != True: #这里要判读得到解是否符合约束条件，如果不符合要重新取一个解
    best_x = get_x(best_x,0,move,t/T)
    jud_n += 1
best_y = target_fun(best_x)
print(f"f = {best_y},x = {best_x},第{jud_n}次循环")

y = best_y

这里定义了一个get_x(x:list,lb,ub,t_move)函数作为获取解的方式，它接收一个当前解x，并对其中的随机一项进行更改，得到新解然后返回。lb和ub是x的变化范围，lb是下限，ub是上限，也就是说x会在lb和ub之间进行变化。不过有时候，我们希望随着温度降低，x的变化范围越来越小，那就可以用t_move影响lb和ub，随着温度的降低t_move可以缩短lb和ub的范围。反之，如果不想x的变化范围受到温度的影响，可以把t_move设为1。

接着重点来了

我们这里需要大概猜一下，x的定义域，也就是x的变化范围。定义域猜的好的话，算法计算的就比较快，如果定义域猜的不好，算法计算的时候就会很慢，很难计算出结果。

看到这两个约束条件，如果x=[0,0,4]是符合条件的，如果x=[0,0,5]就不符合条件了，那就猜个x的定义域是0到5吧，设move等于5，也就是x会在0到5之间随机取一个浮点数作为解。这里只是猜个大概就行。

带入算法进行计算

直接放上完整代码：

import random
import numpy as np
import math
t = 2000
T = 2000
dt = 0.993
eps = 1e-14
n = 0


x = [0,0,0]
y = 100


target_fun = lambda x: x[0]**2 + x[1] ** 2 + x[2]**2 + 8
st_fun1 = lambda x: x[0]**2 - x[1] + x[2]**2 >= 0
st_fun2 = lambda x:  x[0] + x[1]**2 +x[2]**2 <= 20
st_fun3 = lambda x:  math.isclose(-x[0] - x[1]**2 + 2 , 0,rel_tol=0.1, abs_tol=0.001)
st_fun4 = lambda x:  math.isclose(x[1] + 2*x[2]**2 , 3,rel_tol=1e-1, abs_tol=0.01)
st_fun5 = lambda x: x[0] >= 0 and x[1] >= 0 and x[2] >= 0
judgy_st = lambda x: st_fun1(x) and st_fun2(x) and st_fun3(x) and st_fun4(x) and st_fun5(x)


move = 4
t_move = t/T

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

jud_n = 0
best_x = get_x(x,0,move,t/T)
# breakpoint()
while judgy_st(best_x) != True:
    best_x = get_x(best_x,0,move,t/T)
    jud_n += 1
best_y = target_fun(best_x)
print(f"f = {best_y},x = {best_x},第{jud_n}次循环")

y = best_y

while n!= 200 :
    # print(t/T)
    dx = get_x(x,0,move,t/T)
    judgy_n = 0  #这里用一个judgy_n判断内循环的运行次数，如果内循环运行到一定次数还没有计算出解就跳出循环，防止出现死循环
    while judgy_st(dx) != True:
        dx = get_x(dx,0,move,1)
        judgy_n += 1
        if judgy_n == 1000000:#这里设置内循环达到1000000次时跳出循环
            break
    if judgy_n == 1000000:
        continue
    dy = target_fun(dx)

    if dy < y:
        y = dy
        x = dx
        if dy < best_y:
            best_y = dy
            best_x = dx
            print(f"f = {best_y},x = {best_x},第{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)

计算结果如下：

这是视频中给出的解：

可以看到用模拟退火计算的结果与视频的结果还是比较相近的，而且退火算法计算的最小值更小。

总结

这个算法是具有通用性的，只要把目标函数和约束条件改一下，就可以计算类似这样的非线性规划问题，当然也可以解决线性规划问题。