模拟退火算法SAA求解常见问题

模拟退火算法SAA求解常见问题

使用SAA进行数组排序、求解旅行商问题(TSP)、求解背包问题(Knapsack)、多项式回归(Polynomial Regression)。 本文仅仅是为了思考这些问题在SAA算法下如何求解,不考虑求解效率问题。

SAA算法

此算法最小化目标值,如果需要最大化目标值,则需要修改一下公式(1)

Step1: 生成随机解 X, 带入目标函数求解目标值 V. 生成初始温度 T

Step2: 按照一个随机策略,由X生成新的解 X ′ X' X, 并计算新的目标值 V ′ V' V

Step3: 按照公式, 计算出接受新解的概率, P= { 1 , i f V ′ < V e − ( V ′ − V ) / T e l s e \begin{cases}1, \qquad if \quad V'< V \\ e^{-(V'-V) / T} \qquad else\end{cases} {1,ifV<Ve(VV)/Telse   式(1)

Step4: 生成随机数 pr,0 < pr < 1, 且如果 pr < P, 则接收新的解,令 X= X ′ X' X, V= V ′ V' V

Step5: 重复 Step2 - Step4 步骤,N次

Step6: 更新温度, T=T * a, 0< a < 1 或者其他合理的温度下降公式

Step7: 重复 Step2 到 Step6, 直到符合停止条件,例如: 温度小于一个给定的值,或者 V小于一个给定的值, 停止条件可以根据实际情况改变

SAA算法抽象类封装

import abc
import time
import numpy as np
import matplotlib.pyplot as plt

# plt 显示中文设置
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus']=False

class SAA(metaclass=abc.ABCMeta):

    def __init__(self, x=None, max_epoch=float('inf'), t0=0, loop_times=30, rate=0.999):
        self.x = x  # np.ndarray
        if x is not None:
            self.value = self.objective_func(x)
        else:
            self.value = None

        self.best_x = self.copy_x(x)
        self.best_value = self.value

        self.run_data = []  # 存储运行数据
        self.max_epoch = max_epoch  # 最大优化轮次
        self.t = t0  # 初始温度
        self.t0 = t0
        self.epoch = 0
        self.loop_times = loop_times
        self.rate = rate

    def probability(self, delta, t):
        if t <= 0 or delta < 0:
            return 1

        return np.exp(-(delta / t))

    def temperature_down(self, T, rate=0.999, num=0):
        a = self.rate or rate
        #  指数下降温度
        return T * a

    def temperature_init(self):
        if self.t0:
            self.t = self.t0
            return self.t0
        self.t = np.sum(np.array(self.x, 'float64')) ** 2 + 10000
        self.t0 = self.t
        # print("===init_t", self.t)

    @staticmethod
    def copy_x(x):
        return x.copy()

    @abc.abstractmethod
    def update_x(self, x):
        return x

    @abc.abstractmethod
    def objective_func(self, x):
        return 0

    def epoch_callback(self):
        t = self.t
        v = self.value
        bv = self.best_value
        self.run_data.append({"epoch": self.epoch, "temperature": t, "value": v, "best_value": bv})
        return

    def stop(self, min_temprature=0.000001):
        n = self.epoch
        if self.t <= min_temprature or n > self.max_epoch:
            return True
        return False
    
    def show_run_curve(self):
        x = [i['epoch'] for i in self.run_data]
        y =  [i['value'] for i in self.run_data]
        y2 =  [i['best_value'] for i in self.run_data]
        plt.plot(x, y, label='value')
        plt.plot(x, y2, label='best_value')
        plt.title('优化曲线')
        plt.legend()
        plt.show()

    def run(self, clear_run_data=True):
        self.temperature_init()
        self.epoch = 0
        loop_times = self.loop_times
        
        print(f"初始温度:{self.t}, loop_times: {loop_times}")
        if clear_run_data:
            self.run_data.clear()
            
        if self.value is None:
            self.value = self.objective_func(self.x)
            self.best_value = self.value

        best_value = self.best_value
        while True:
            for _ in range(loop_times):
                x = self.x
                # old_value = self.objective_func(x)
                old_value = self.value
                new_x = self.update_x(x)
                new_value = self.objective_func(new_x)
                
                delta = new_value - old_value
                p = self.probability(delta, self.t)
                pr = np.random.random()
    
                if best_value > new_value:  # 沿着value 下降的方向优化
                    self.best_x = self.copy_x(new_x)
                    self.best_value = new_value
                    best_value = new_value
    
                if pr < p:  # 按概率更新解
                    self.x = new_x
                    self.value = new_value
                    # print(f"==={n}概率接受解p: {p} T: {self.t} δ:{delta} ({pr},{p})")

            print(f"==epoch: {self.epoch} t:{self.t} {old_value} -> {new_value}", end='\r')
            # 回调
            self.epoch_callback()
            # 温度下降
            self.t = self.temperature_down(self.t, num=self.epoch)

            # 停机判断
            if self.stop():
                print(f"===epoch: {self.epoch} 退出:T:{self.t}, P:{p}, {old_value}")
                break

            self.epoch += 1

        self.x = x
    
    def run_with_t0_update(self, times=5, show_curve=True):
        for i in range(times):
            self.t0 *= 1 / (10 ** i)
            self.run()
            self.show_run_curve()
            # show_path(tsp102.best_x, points_xy2)

排序

定义从小到大的排序目标函数为:左边大于右边的个数. 优化目标函数达到最小值.

import numpy as np


class SORT(SAA):
    def update_x(self, x):
        x = x.copy()
        i0 = np.random.randint(0, x.shape[-1])
        i1 = np.random.randint(0, x.shape[-1])
        while x.shape[-1] > 2 and i0 == i1:
            i1 = np.random.randint(0, x.shape[-1])
            
        x[i0], x[i1] = x[i1], x[i0]
        return x
    
    def objective_func(self, x):
        n = 0
        for i in range(x.shape[-1] - 1):
            n += int(x[i+1] < x[i])
        
        time.sleep(0.00000001)
        return n

    def stop(self, min_temprature=0.000001):
        # 排序问题的最优值是可知的, 可以提前退出计算
        if self.best_value <= 0:
            return True
        
        n = self.epoch
        if self.t <= min_temprature or n > self.max_epoch:
            return True
        return False
s = SORT(x=np.array([0, 1, 2, 3, 1,3233, 2, 29, 10, 889, 9]), t0=100)
# s = SORT(x=np.array([600000000, 5, 1, 2, 3, -4]))
# s = SORT(x=np.array([60000, 1, 2]))
s.run()
初始温度:100, loop_times: 30
===epoch: 5672 退出:T:0.3427793934809477, P:0.0029415406738396072, 0
s.show_run_curve()

在这里插入图片描述

print(list(s.x))
[0, 1, 1, 2, 2, 3, 9, 10, 29, 889, 3233]
print(list(s.best_x))
[0, 1, 1, 2, 2, 3, 9, 10, 29, 889, 3233]

TSP问题求解

TSP问题,即旅行商问题(Traveling Salesman Problem),是数学领域中的一个著名问题。
它是一个组合优化问题,旨在找到访问给定城市集合的最短可能路径,同时满足每个城市只能访问一次,并且最终要回到起点。
TSP问题可以用以下方式定义:

问题定义:给定一个包含n个城市的集合,找到访问每个城市一次并最终回到起点的最短路径。

目标:使得所走的路程最短。

限制:每个城市只能访问一次。

TSP问题是一个NP-完全问题(NPC),这意味着随着城市数量的增加,问题的计算复杂性会显著增加,使得在合理时间内找到精确解变得非常困难。
此外,TSP问题还有多个变体,包括对称TSP和非对称TSP,以及带有时间窗口和资源限制的TSP问题。

import numpy as np

class TSP(SAA):

    def __init__(self, matrix, points=None, *args, **kwargs):
        self.matrix = matrix
        self.points = points
        super().__init__(*args, **kwargs)
        
        
    def update_x(self, x):
        x = x.copy()
        i0 = np.random.randint(1, x.shape[-1])
        i1 = np.random.randint(1, x.shape[-1])
        while x.shape[-1] > 3 and i0 == i1:
            i1 = np.random.randint(1, x.shape[-1])
            
        x[i0], x[i1] = x[i1], x[i0]
        return x
    
    def objective_func(self, x):
        n = self.matrix[x[0], x[-1]]
        for i in range(x.shape[-1] - 1):
            n += self.matrix[x[i], x[i+1]]
        return n

    @staticmethod
    def x_to_chars(x, points=''):
        # points = self.points or points
        return [points[i] for i in x]
    
    @staticmethod
    def generate_matrix(points='', distance_matrix=True):
        d = len(points)
        mat = np.random.random((d, d))
        if distance_matrix:
            # 使用numpy.triu生成上三角矩阵
            mat += np.triu(mat, k=1)
            mat += mat.T
            # 填充对角线上的值
            np.fill_diagonal(mat, 0)
            return mat
        return mat

    @staticmethod
    def generate_matrix_by_points(n=6, width=10):
        x = np.random.random(n) *  width
        y = np.random.random(n) *  width
        plt.scatter(x, y)
        plt.show()
        points = list(zip(x, y))
        mat = np.zeros((n, n), dtype='float32')
        for i in range(n):
            for j in range(i + 1, n):
                d = ((points[i][0]  - points[j][0]) ** 2 + (points[i][1]  - points[j][1]) ** 2) ** 0.5
                mat[i][j] = d
                mat[j][i] = d
        return mat, points

    def show_run_curve(self):
        plt.rcParams['font.sans-serif'] = ['SimHei']
        x = [i['epoch'] for i in self.run_data]
        y =  [i['value'] for i in self.run_data]
        y2 =  [i['best_value'] for i in self.run_data]
        plt.plot(x, y, label='value')
        plt.plot(x, y2, label='best_value')
        plt.title('TSP优化曲线')
        plt.legend()
        plt.show()

    def run_with_t0_update(self, times=5, show_curve=True, points_xy=None):
        for i in range(times):
            self.t0 *= 1 / (10 ** i)
            self.run()
            self.show_run_curve()
            if points_xy:
                show_path(self.best_x, points_xy)
            
points = list('abcdefghjklmn')
# m = {(i,j): np.random.random() *  100 for i, _ in enumerate(points) for j, _ in enumerate(points)}
# print(m)
np.random.seed(200)
m = TSP.generate_matrix(points=points)
# print(m)
m, points_xy = TSP.generate_matrix_by_points(n = len(points))

在这里插入图片描述

tsp = TSP(m, points=points, x=np.array(range(len(points))))
tsp.run()
初始温度:16084.0, loop_times: 30
===epoch: 23489 退出:T:9.99338209184425e-07, P:0.0, 34.03287887573242417968750164
tsp.x_to_chars(tsp.x, points=tsp.points), tsp.x_to_chars(tsp.best_x, points=tsp.points)
(['a', 'b', 'e', 'j', 'g', 'm', 'n', 'h', 'c', 'k', 'l', 'd', 'f'],
 ['a', 'd', 'l', 'k', 'c', 'h', 'n', 'm', 'g', 'j', 'e', 'b', 'f'])

优化曲线

tsp.show_run_curve()

在这里插入图片描述

def show_path(path, points):
    x = [i[0] for i in points]
    y = [i[1] for i in points]
    plt.scatter(x, y)
    path_x = [points[i][0] for i in path]
    path_y = [points[i][1] for i in path]
    # 绘制箭头
    s = 0.18
    for i in range(0, len(path) - 1):
        x_start, y_start = path_x[i], path_y[i]
        x_end, y_end = path_x[i + 1], path_y[i + 1]
        
        a = x_end - x_start
        b = y_end - y_start
        if a == 0 and b == 0:
            continue
            
        l = 1 - s / ((a ** 2 + b ** 2) ** 0.5)
        a *= l
        b *= l
        plt.arrow(x_start, y_start, 
                  a, b, 
                  head_width=0.2, head_length=s, fc='red', ec='red')
    
    # plt.plot(path_x, path_y, 'r--')
    plt.show()

show_path(tsp.best_x, points_xy)

在这里插入图片描述

2-opt方法搜索

class TSPByKopt(TSP):
    def update_x(self, x):
       return self.opt2(x)

    @staticmethod
    def opt2(x):
        x = x.copy()
        i0 = np.random.randint(1, x.shape[-1])
        i1 = np.random.randint(1, x.shape[-1])
        while x.shape[-1] > 3 and i0 == i1:
            i1 = np.random.randint(1, x.shape[-1])
        
        while i0 < i1:
            x[i0], x[i1] = x[i1], x[i0]
            i0 += 1
            i1 -= 1
        return x
tsp2 = TSPByKopt(m, points=points,x=np.array(range(len(points))))
tsp2.run()
初始温度:16084.0
===epoch: 23489 退出:T:9.99338209184425e-07, P:1.0, 33.72549819946289498199462896
tsp2.x_to_chars(tsp2.x, points=tsp2.points), tsp2.x_to_chars(tsp2.best_x, points=tsp2.points)
(['a', 'd', 'l', 'k', 'c', 'h', 'n', 'm', 'g', 'j', 'e', 'b', 'f'],
 ['a', 'd', 'l', 'k', 'c', 'h', 'n', 'm', 'g', 'j', 'e', 'b', 'f'])
tsp2.show_run_curve()

在这里插入图片描述

show_path(tsp2.best_x, points_xy)

在这里插入图片描述

points = list(range(100))
m2, points_xy2 = TSP.generate_matrix_by_points(n = len(points))
tsp101 = TSP(m2, points=points,x=np.array(range(len(points))))

在这里插入图片描述

两点对换与2-opt方法比较

tsp101.run()
初始温度:24512500.0
===epoch: 30814 退出:T:9.997865508578958e-07, P:0.0, 145.0625610351562540283203125
tsp101.show_run_curve()

在这里插入图片描述

show_path(tsp101.best_x, points_xy2)

在这里插入图片描述

points = list(range(100))
# m2, points_xy2 = TSPByKopt.generate_matrix_by_points(n = len(points))
tsp102 = TSPByKopt(m2, points=points,x=np.array(range(len(points))))
tsp102.run()
初始温度:24512500.0
===epoch: 30814 退出:T:9.997865508578958e-07, P:1.0, 105.5768737792968873779296888
tsp102.show_run_curve()
show_path(tsp102.best_x, points_xy2)

在这里插入图片描述

在这里插入图片描述


tsp102.best_value, tsp102.best_value
(105.57687, 105.57687)
tsp102.run_with_t0_update()
初始温度:0.00245125
===epoch: 7800 退出:T:9.99450426296383e-07, P:0.0, 93.37339019775395860290527343

在这里插入图片描述

在这里插入图片描述

初始温度:0.000245125
===epoch: 5499 退出:T:9.990169333340238e-07, P:1.0, 92.5606842041015684204101566

在这里插入图片描述

在这里插入图片描述

初始温度:2.45125e-06
===epoch: 896 退出:T:9.991496610472398e-07, P:1.0, 92.5606842041015642041015675

求解背包问题

我们有n种物品,物品j的重量为wj,价格为pj

goods = 10 # 总物品数
np.random.seed(10)
# 需要求解的变量, 初始变量
selects = np.random.randint(0, 2, size=10)
# 每个物品可以被选择的次数
select_times = np.random.randint(1, 5, size=10)
weights = np.random.random(goods) # 物品重量
values = np.random.randint(1, 20, size=10) # 物品价值
total = 10 # 背包总容量
print(weights)
print(values)
print(selects)
print(f"select_times {select_times}")
[0.68535982 0.95339335 0.00394827 0.51219226 0.81262096 0.61252607
 0.72175532 0.29187607 0.91777412 0.71457578]
[14 14 13  2  5 19 14 12 11 10]
[1 1 0 1 0 1 1 0 1 1]
select_times [1 2 2 3 1 2 1 3 1 3]

class Knapsack(SAA):
    def __init__(self, weights, values, select_times, total, x=None, *args, **kwargs):
        self.weights = weights
        self.values=values
        self.select_times = select_times
        self.total = total
        if x is not None:
            self.constraint(x)
        super().__init__(x=x, *args, **kwargs)

    # def probability(self, delta, t):
    #     if t <= 0 or delta < 0:
    #         return 1

    #     return np.exp(-(delta / t))
        
    def update_x(self, x):
        x = x.copy()
        self.constraint(x)
        
        i = np.random.randint(0, len(x))
        if np.random.random() > 0.5:
            times = np.random.randint(0, self.select_times[i] + 1)
            x[i] = 0
            x[i] = max(min(int(self.total - x.sum()), times), 0)
        else:
            d = (np.random.random() > 0.5) * 2 - 1
            a = min(self.select_times[i], x[i] + d)
            x[i] = 0
            x[i] = max(min(int(self.total - x.sum()), a), 0)
        
        if x[i] > self.select_times[i]:
            print("解异常")
            
        self.constraint(x)
        return x

    def constraint(self, x):
        """
        背包问题约束条件
        """
        for i in range(len(x)):
            while (x[i] > 0 and x.sum() > self.total) or (x[i] > self.select_times[i]):
                x[i] = max(0, x[i] - 1)
            
            if x[i] > self.select_times[i]:
                raise Exception("约束异常")
                
    def objective_func(self, x):
        v = 0
        for i in range(len(x)):
            v += x[i] * self.values[i]

        # 转换成求最小值
        m = self.values.sum() * self.select_times.sum()
        return m - v
        # return -v
        
knapsack = Knapsack(weights, values, select_times, total, x=selects)
knapsack.run()
knapsack.show_run_curve()
初始温度:10049.0
===epoch: 23019 退出:T:9.992215880401829e-07, P:1.0, 2032

在这里插入图片描述

print(knapsack.best_x,knapsack.best_value)
print(knapsack.select_times)
[1 2 2 0 0 2 1 2 0 0] 2022
[1 2 2 3 1 2 1 3 1 3]

最短路径

s_points = list(range(50))
# shortest path problem
class SPP(TSP):
    """
    邻接矩阵为无穷大表示不连通
    """
    def update_x(self, x):
        # assert x.shape[-1] > 3, "点数量过少"
        x = x.copy()
        self.constraint(x)

        if np.random.random() > 0.6:
            i0 = np.random.randint(1, x.shape[-1] - 1) # 起止点为0 和 -1, 不交换
            i1 = np.random.randint(1, x.shape[-1] - 1)
            while x.shape[-1] > 4 and i0 == i1:
                i1 = np.random.randint(1, x.shape[-1] - 1)

            if np.random.random() > 0.5:
                x[i0], x[i1] = x[i1], x[i0]
            else:
                # 直接取消中间的路径, 如果不加此方法,则很难找到最优解
                while i0 < i1:
                    x[i0 + 1] = x[i0]
                    i0 += 1
                    
        else:
            i0 = np.random.randint(1, x.shape[-1] - 1)
            s = np.random.randint(0, self.matrix.shape[-1] - 1)
            x[i0] = s
            
        self.constraint(x)
        return x
    
    def objective_func(self, x):
        n = 0
        inf = float('inf')
        # print(self.matrix)
        for i in range(x.shape[-1] - 1):
            if self.matrix[x[i], x[i+1]] == inf:
                n += np.abs(self.matrix).sum()
                
            n += self.matrix[x[i], x[i+1]]
        return n
    
    def constraint(self, x):
        """
        约束条件
        """
        for i in range(x.shape[-1]- 2):
            if self.matrix[x[i], x[i + 1]] == float('inf'):
                x[i + 1] = x[i].copy()
np.random.seed(6)
m, p_list = SPP.generate_matrix_by_points(n=len(s_points))
x = np.array(range(len(s_points)))
np.random.shuffle(x)
max_xy = np.argwhere(m == m.max())[0]
print(max_xy)
x[0], x[-1] = max_xy
print(x)
# 起止点设置不连通
m[x[0], x[-1]] = float('inf')
shortest = SPP(matrix=m, points=s_points, x=x)

在这里插入图片描述

[31 32]
[31 10 24 38 28 33 11  4 22 13 46 17 12 37 40 34 21 36  7 27  9 26  1 18
 48  6 47  0 25 30 29 20 14  2 43 41 16 31 45 23 32 44 49 42 35 39 15  8
  3 32]
# shortest.t0 = 10 ** 3
shortest.run_with_t0_update(times=2, points_xy=p_list)
# shortest.show_run_curve()
初始温度:1607696.0
===epoch: 28091 退出:T:9.99767870003539e-07, P:0.0, 11.777696013450623181731224065

在这里插入图片描述

在这里插入图片描述

初始温度:160769.6
===epoch: 25790 退出:T:9.993342393559143e-07, P:1.0, 11.77769601345062396013450623

退火算法查找最短路径与Dijkstra算法比较

def dijkstra(matrix, i0=0):
    n = matrix.shape[0]
    find = [-1 for i in range(n)]
    distance = [float("inf") for i in range(n)]
    pre = [-1 for i in range(n)]
    
    distance[i0] = 0
    find[i0] = 1
    for j in range(n):
        if j != i0:
            distance[j] = matrix[i0, j]
            pre[j] = i0

    # print(distance)
    
    i = i0
    while True:
        # 遍历每一个未确认的点,尝试更新未确认点的距离,即:通过未确认的点,则距离变短
        for j in range(n):
            if find[j] > 0 or i == j:
                continue
                
            if matrix[i,j] >= float("inf"):
                continue
            
            new_d = matrix[i,j] + distance[i]
            if new_d < distance[j]: # 当从i经过j点时,距离可以更短,则更新j的最小距离
                distance[j] = new_d
                pre[j] = i

        # 在未确认的点中找到最小距离的点
        m = float("inf")
        i = None
        for j in range(n):
            if find[j] > 0:
                continue
            if distance[j] < m:
                i = j
                m = distance[j]
                
        if i is None:
            break
        
        # print(f"next: {i}")
        find[i] = 1

    print('=' * 10)
    print(find)
    print(distance)
    print(pre)
    return distance, pre

def to_path(s, last=0):
    p = [last]
    while s[last] != -1:
        last = s[last]
        p.insert(0, last)
    
    return p
distance, pre_points = dijkstra(m, i0=x[0])
dijkstra_path = to_path(pre_points, last=x[-1])
==========
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
[7.304405, 7.0522614, 2.2776055, 10.450214, 9.338489, 8.000453, 4.8732204, 6.442265, 10.385517, 6.835383, 10.106222, 9.702004, 9.271525, 4.3813324, 6.7429028, 4.347516, 4.0250196, 9.10099, 7.93373, 6.641005, 9.448183, 7.1655736, 1.9282913, 2.9891455, 9.282758, 5.8116393, 8.759328, 6.4195356, 6.034359, 9.861004, 3.476944, 0, 11.777696, 3.2385242, 6.152703, 4.0818205, 6.212417, 10.283808, 9.299705, 8.050401, 2.5899367, 7.2155046, 6.103101, 4.4901824, 5.4372907, 6.782229, 8.695412, 10.481743, 7.812679, 3.9552972]
[31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, -1, 37, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31, 31]
print(dijkstra_path)
[31, 37, 32]
show_path(dijkstra_path, p_list)

在这里插入图片描述

多项式回归(Polynomial Regression)

def polynomial_func(x, params):
    """
    params type: [(float, float)]
    """
    v = 0
    for a, b in params:
        v += a * np.power(x, b)
    
    return v
## 生成一个随机的多项式函数,测试代码
pa = [(np.random.random(),  i) for i in range(10)]
pa = np.array(pa)
pa, pa.shape
(array([[0.17746591, 0.        ],
        [0.89390912, 1.        ],
        [0.32112076, 2.        ],
        [0.12548015, 3.        ],
        [0.34522942, 4.        ],
        [0.69496884, 5.        ],
        [0.02675415, 6.        ],
        [0.5419919 , 7.        ],
        [0.29978091, 8.        ],
        [0.66899116, 9.        ]]),
 (10, 2))
x = np.arange(-2, 2, 0.001)
y = np.apply_along_axis(lambda a: polynomial_func(a, pa), 0, x)
plt.plot(x, y)
plt.show()

在这里插入图片描述

拟合二次函数

np.random.seed(10)
train_x = np.arange(-2, 2,  0.01)
train_y = -np.power(train_x, 2) + np.random.random(len(train_x)) * 0.3

plt.scatter(train_x, train_y)
plt.show()

在这里插入图片描述

class PolynomialRegression(SAA):
    def __init__(self, train_x, train_y, *args, **kwargs):
        self.train_x = train_x
        self.train_y = train_y
        
        super().__init__(*args, **kwargs)
        
    def update_x(self, x):
        i = np.random.randint(0, x.shape[0])
        if np.random.random() < 0.01:
            x[i, 0] = np.random.random() * 2 - 1
        else: 
            a = (np.random.random() * 2 - 1) * 0.1
            x[i, 0] *= a
        return x
        
    def objective_func(self, x):
        pre_y = np.apply_along_axis(lambda a: polynomial_func(a, x), 0, self.train_x)
        d = pre_y - self.train_y
        
        return (d ** 2).sum() / self.train_y.shape[0]
        
    def predict(self, x, params=None):
        if params is None:
            params = self.best_x
        
        return np.apply_along_axis(lambda a: polynomial_func(a, params), 0, x)
pol = PolynomialRegression(train_x, train_y, x=pa)
pol.run_with_t0_update(times=5)
初始温度:12410.3870039805, loop_times: 30
===epoch: 23230 退出:T:9.991767101361876e-07, P:0.0, 0.0171797724867069288576094586636

在这里插入图片描述

初始温度:1241.0387003980502, loop_times: 30
===epoch: 20928 退出:T:9.99743078972076e-07, P:0.0, 0.0169463164575696669857609458566

在这里插入图片描述

初始温度:12.410387003980503, loop_times: 30
===epoch: 16325 退出:T:9.99875903159779e-07, P:0.0, 0.0170044256091905779839097938762

在这里插入图片描述

初始温度:0.012410387003980504, loop_times: 30
===epoch: 9421 退出:T:9.99575009871613e-07, P:0.0, 0.01694526861586353778185395333466

在这里插入图片描述

初始温度:1.2410387003980505e-06, loop_times: 30
===epoch: 215 退出:T:9.998406312320972e-07, P:0.0, 0.01694526861586353757609458666

在这里插入图片描述

pre_y = pol.predict(train_x)
np.round(pol.best_x, 2)
array([[ 0.  ,  0.  ],
       [-0.  ,  1.  ],
       [-0.94,  2.  ],
       [ 0.  ,  3.  ],
       [ 0.  ,  4.  ],
       [-0.  ,  5.  ],
       [-0.  ,  6.  ],
       [-0.  ,  7.  ],
       [-0.  ,  8.  ],
       [ 0.  ,  9.  ]])
plt.scatter(train_x, train_y, c='g')
plt.scatter(train_x, pre_y, c='r')
plt.show()

在这里插入图片描述

拟合三次函数

np.random.seed(10)
train_x = np.arange(-2, 2,  0.01)
train_y = np.power(train_x, 3) + np.random.random(len(train_x))

plt.scatter(train_x, train_y)
plt.show()

在这里插入图片描述

pol = PolynomialRegression(train_x, train_y, x=pa)
pol.run()
初始温度:12024.999999982396, loop_times: 30
===epoch: 23198 退出:T:9.996464766364758e-07, P:0.0, 0.20693877934235541598321377945
np.round(pol.best_x, 2)
array([[ 0.46,  0.  ],
       [-0.  ,  1.  ],
       [-0.  ,  2.  ],
       [ 0.88,  3.  ],
       [-0.  ,  4.  ],
       [ 0.  ,  5.  ],
       [-0.  ,  6.  ],
       [-0.  ,  7.  ],
       [ 0.  ,  8.  ],
       [ 0.  ,  9.  ]])
pol.show_run_curve()

在这里插入图片描述

pre_y = pol.predict(train_x)
plt.scatter(train_x, train_y, c='g')
plt.scatter(train_x, pre_y, c='r')
plt.show()

在这里插入图片描述

  • 8
    点赞
  • 22
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值