对于数学建模里用到的多目标优化算法进行探究
使用上次建模选拔赛里做的匆忙且不完整的调参作为优化问题
选择对下面的方程进行拟合确定参数
选择 [1:100] 的100个散点作为拟合区间

def realFuction(x):
    return -0.02*x**2 +2*x+5*math.sin(x)
def myFuc(a1,a2,a3,x):
    return  a1*x**2 + a2*x + a3*math.sin(x)

模拟退火

基于参考链接的代码
但是源代码我认为有一些问题：第57行原代码为rp，y和yNew相近的时候我认为逼近优解，所以应该用小概率接受较差解，所以改了。感觉效果也有提升

相对于源代码，我加上了保存最优解的代码，在后面的结果展示里也可以看到对应的结果，防止在最后的遍历里小概率的选择了较差解，导致结果极差

模拟退火参数选择：

1. 对于T0决定了一开始参数的变化幅度（参数变化乘了T），所以不易选择过大（容易震荡）或过小
2. T0和Tmin又共同决定了温度变化的幅度，影响了参数的变化幅度和退火过程长度
3. T = T0 / math.log2(1 + t)决定温度变化速度，不加log更快，但是结果不一定更好
4. 对于k即每一个温度下对参数的尝试调整次数 ```python from future import division import numpy as np import matplotlib.pyplot as plt import math from Draw import Draw def realFuction(x): return -0.02x**2 +2x+5*math.sin(x)

def myFuc(a1,a2,a3,x): return a1x**2 + a2x + a3*math.sin(x)

def aimFunction(para,x): res=0 for i in x: res += math.sqrt((myFuc(para[0],para[1],para[2],i)-realFuction(i))**2) return res/len(x)

x = [i for i in range(100)] y = [0 for i in range(100)] for i in range(100): y[i] = realFuction(x[i])

T0 = 200 # initiate temperature\ T=T0 Tmin = 15 # minimum value of terperature

para=[np.random.uniform(low=-1, high=5) for _ in range(3)]

k = 100 # times of internal circulation temp_y = 0 # initiate result t = 0 # time best=[] best_fit=0 while T >= Tmin: for i in range(k):

    # calculate y
    temp_y = aimFunction(para,x)
    # generate a new x in the neighboorhood of x by transform function
    paraNew = [i+np.random.uniform(low=-0.0005,high=0.0005)*T for i in para]
    yNew = aimFunction(paraNew,x)
    if len(best)==0:
        best_fit=yNew
        best=paraNew
    else:
        if yNew<best_fit:
            best=paraNew
            best_fit=yNew
    if yNew - temp_y < 0:
        para = paraNew
    else:
        # metropolis principle
        p = math.exp(-(yNew - temp_y) / T)
        r = np.random.uniform(low=0, high=1)
        if r > p:
            para = paraNew
print(para,aimFunction(para,x),best_fit)
t += 1
T = T0 / math.log2(1 + t)
print("cur T: ",T)

result=[] bestresult=[] for i in x: result.append(myFuc(para[0],para[1],para[2],i)) bestresult.append(myFuc(best[0],best[1],best[2],i)) print(para) print(best) Draw([y,result,bestresult],[“real”,”predict”,”best”],x,”Real and predict”,Colors=None)

（效果比较好的图忘记截图参数结果了）<br />![Figure_1.png](https://cdn.nlark.com/yuque/0/2020/png/296244/1598551674668-4babb2e7-5f61-4314-b47a-9444d1112fd2.png#align=left&display=inline&height=700&margin=%5Bobject%20Object%5D&name=Figure_1.png&originHeight=700&originWidth=1200&size=77340&status=done&style=none&width=1200)<br />（重新跑了，调了一些参数的结果，Best的拟合效果巨好，但是perdict飘了，可以看截图里的最后一次退火计算，不知道为啥第一个参数突变了，但是Best效果实在太好了）<br />![image.png](https://cdn.nlark.com/yuque/0/2020/png/296244/1598551698610-1691f3fd-a0e5-4724-b710-7542441004f5.png#align=left&display=inline&height=408&margin=%5Bobject%20Object%5D&name=image.png&originHeight=817&originWidth=1341&size=114057&status=done&style=none&width=670.5)<br />![image.png](https://cdn.nlark.com/yuque/0/2020/png/296244/1598551711489-74d375c9-7bd8-4699-afb9-dbfe2faba449.png#align=left&display=inline&height=142&margin=%5Bobject%20Object%5D&name=image.png&originHeight=284&originWidth=1090&size=38053&status=done&style=none&width=545)
<a name="j8Uj0"></a>
# 遗传算法
基于参考资料做了一下修改<br />构造了问题实例和损失函数，同时修改了突变的操作<br />（原来的是在范围里直接随机一个，实在感觉有点随便，对于解空间比较大的参数可能不适合？）<br />最后修改了Fitness的计算<br />原先是fit值向大的方向优化，这里是需要缩小gap<br />第一次改为负号，不可以，会影响到选择的时候（基于Fitness占比选择遗传，负数导致错误）<br />第二次改为倒数
```python
import random
import math
import numpy as np
from operator import itemgetter
def realFuction(x):
    return -0.02*x**2 +2*x+5*math.sin(x)
def myFuc(a1,a2,a3,x):
    return  a1*x**2 + a2*x + a3*math.sin(x)
def aimFunction(para,x):
    res=0
    for i in x:
        res += math.sqrt((myFuc(para[0],para[1],para[2],i)-realFuction(i))**2)
    return 1/(res/len(x))
x = [i  for i in range(100)]
y = [0 for i in range(100)]
for i in range(100):
    y[i] = realFuction(x[i])
class Gene:
    """
    This is a class to represent individual(Gene) in GA algorithom
    each object of this class have two attribute: data, size
    """
    def __init__(self, **data):
        self.__dict__.update(data)
        self.size = len(data['data'])  # length of gene
class GA:
    """
    This is a class of GA algorithm.
    """
    def __init__(self, parameter):
        """
        Initialize the pop of GA algorithom and evaluate the pop by computing its' fitness value.
        The data structure of pop is composed of several individuals which has the form like that:
        {'Gene':a object of class Gene, 'fitness': 1.02(for example)}
        Representation of Gene is a list: [b s0 u0 sita0 s1 u1 sita1 s2 u2 sita2]
        """
        # parameter = [CXPB, MUTPB, NGEN, popsize, low, up]
        self.parameter = parameter
        low = self.parameter[4]
        up = self.parameter[5]
        self.bound = []
        self.bound.append(low)
        self.bound.append(up)
        pop = []
        for i in range(self.parameter[3]):
            geneinfo = []
            for pos in range(len(low)):
                geneinfo.append(random.randint(self.bound[0][pos], self.bound[1][pos]))  # initialise popluation
            fitness = self.evaluate(geneinfo)  # evaluate each chromosome
            pop.append({'Gene': Gene(data=geneinfo), 'fitness': fitness})  # store the chromosome and its fitness
        self.pop = pop
        self.bestindividual = self.selectBest(self.pop)  # store the best chromosome in the population
    def evaluate(self, geneinfo):
        """
        fitness function
        """
        y = aimFunction(geneinfo,x)
        return y
    def selectBest(self, pop):
        """
        select the best individual from pop
        """
        s_inds = sorted(pop, key=itemgetter("fitness"), reverse=True)  # from large to small, return a pop
        return s_inds[0]
    def selection(self, individuals, k):
        """
        select some good individuals from pop, note that good individuals have greater probability to be choosen
        for example: a fitness list like that:[5, 4, 3, 2, 1], sum is 15,
        [-----|----|---|--|-]
        012345|6789|101112|1314|15
        we randomly choose a value in [0, 15],
        it belongs to first scale with greatest probability
        """
        s_inds = sorted(individuals, key=itemgetter("fitness"),
                        reverse=True)  # sort the pop by the reference of fitness
        sum_fits = sum(ind['fitness'] for ind in individuals)  # sum up the fitness of the whole pop
        chosen = []
        for i in range(k):
            u = random.random() * sum_fits  # randomly produce a num in the range of [0, sum_fits], as threshold
            sum_ = 0
            for ind in s_inds:
                sum_ += ind['fitness']  # sum up the fitness
                if sum_ >= u:
                    # when the sum of fitness is bigger than u, choose the one, which means u is in the range of
                    # [sum(1,2,...,n-1),sum(1,2,...,n)] and is time to choose the one ,namely n-th individual in the pop
                    chosen.append(ind)
                    break
        # from small to large, due to list.pop() method get the last element
        chosen = sorted(chosen, key=itemgetter("fitness"), reverse=False)
        return chosen
    def crossoperate(self, offspring):
        """
        cross operation
        here we use two points crossoperate
        for example: gene1: [5, 2, 4, 7], gene2: [3, 6, 9, 2], if pos1=1, pos2=2
        5 | 2 | 4  7
        3 | 6 | 9  2
        =
        3 | 2 | 9  2
        5 | 6 | 4  7
        """
        dim = len(offspring[0]['Gene'].data)
        geninfo1 = offspring[0]['Gene'].data  # Gene's data of first offspring chosen from the selected pop
        geninfo2 = offspring[1]['Gene'].data  # Gene's data of second offspring chosen from the selected pop
        if dim == 1:
            pos1 = 1
            pos2 = 1
        else:
            pos1 = random.randrange(1, dim)  # select a position in the range from 0 to dim-1,
            pos2 = random.randrange(1, dim)
        newoff1 = Gene(data=[])  # offspring1 produced by cross operation
        newoff2 = Gene(data=[])  # offspring2 produced by cross operation
        temp1 = []
        temp2 = []
        for i in range(dim):
            if min(pos1, pos2) <= i < max(pos1, pos2):
                temp2.append(geninfo2[i])
                temp1.append(geninfo1[i])
            else:
                temp2.append(geninfo1[i])
                temp1.append(geninfo2[i])
        newoff1.data = temp1
        newoff2.data = temp2
        return newoff1, newoff2
    def mutation(self, crossoff, bound):
        """
        mutation operation
        """
        '''
        dim = len(crossoff.data)
        if dim == 1:
            pos = 0
        else:
            pos = random.randrange(0, dim)  # chose a position in crossoff to perform mutation.
        #
        crossoff.data[pos] = random.uniform(bound[0][pos], bound[1][pos])
        '''
        crossoff.data = [i + np.random.uniform(low=-0.001, high=0.001) for i in crossoff.data]
        return crossoff
    def GA_main(self):
        """
        main frame work of GA
        """
        popsize = self.parameter[3]
        print("Start of evolution")
        # Begin the evolution
        for g in range(NGEN):
            print("------ Generation {} ------".format(g))
            # Apply selection based on their converted fitness
            selectpop = self.selection(self.pop, popsize)
            nextoff = []
            while len(nextoff) != popsize:
                # Apply crossover and mutation on the offspring
                # Select two individuals
                offspring = [selectpop.pop() for _ in range(2)]
                if random.random() < CXPB:  # cross two individuals with probability CXPB
                    crossoff1, crossoff2 = self.crossoperate(offspring)
                    if random.random() < MUTPB:  # mutate an individual with probability MUTPB
                        muteoff1 = self.mutation(crossoff1, self.bound)
                        muteoff2 = self.mutation(crossoff2, self.bound)
                        fit_muteoff1 = self.evaluate(muteoff1.data)  # Evaluate the individuals
                        fit_muteoff2 = self.evaluate(muteoff2.data)  # Evaluate the individuals
                        nextoff.append({'Gene': muteoff1, 'fitness': fit_muteoff1})
                        nextoff.append({'Gene': muteoff2, 'fitness': fit_muteoff2})
                    else:
                        fit_crossoff1 = self.evaluate(crossoff1.data)  # Evaluate the individuals
                        fit_crossoff2 = self.evaluate(crossoff2.data)
                        nextoff.append({'Gene': crossoff1, 'fitness': fit_crossoff1})
                        nextoff.append({'Gene': crossoff2, 'fitness': fit_crossoff2})
                else:
                    nextoff.extend(offspring)
            # The population is entirely replaced by the offspring
            self.pop = nextoff
            # Gather all the fitnesses in one list and print the stats
            fits = [ind['fitness'] for ind in self.pop]
            best_ind = self.selectBest(self.pop)
            if best_ind['fitness'] > self.bestindividual['fitness']:
                self.bestindividual = best_ind
            print("Best individual found is {}, {}".format(self.bestindividual['Gene'].data,self.bestindividual['fitness']))
            print("Max fitness of current pop: {}".format(max(fits)))
        print("------ End of (successful) evolution ------")
from Draw import Draw
if __name__ == "__main__":
    CXPB, MUTPB, NGEN, popsize = 0.7, 0.5, 1000, 1000#control parameters
    up = [5, 5, 5]  # upper range for variables
    low = [-2, -2, -2]  # lower range for variables
    parameter = [CXPB, MUTPB, NGEN, popsize, low, up]
    run = GA(parameter)
    run.GA_main()
    print(run.bestindividual['Gene'].data)
    result = []
    for i in x:
        result.append(myFuc(run.bestindividual['Gene'].data[0], run.bestindividual['Gene'].data[1], run.bestindividual['Gene'].data[2], i))
    Draw([y, result], ["real", "predict"], x, "Real and predict", Colors=None)

蒙特卡洛算法

还有一种算法就是使用蒙特卡洛进行优化
优化的方法很简单，就是直接随机的生成参数，然后进行尝试，如果优化则保留
简单的基于模拟退火的算法框架实现如下

from __future__ import division
import numpy as np
import matplotlib.pyplot as plt
import math
from Draw import Draw
def realFuction(x):
    return -0.02*x**2 +2*x+5*math.sin(x)
def myFuc(a1,a2,a3,x):
    return  a1*x**2 + a2*x + a3*math.sin(x)
def aimFunction(para,x):
    res=0
    for i in x:
        res += math.sqrt((myFuc(para[0],para[1],para[2],i)-realFuction(i))**2)
    return res/len(x)
x = [i  for i in range(100)]
y = [0 for i in range(100)]
for i in range(100):
    y[i] = realFuction(x[i])
Need=[np.random.uniform(low=-1, high=5) for _ in range(3)]
cur_best=aimFunction(Need,x)
update_time=0
for i in range(100000):
    para=[np.random.uniform(low=-1, high=5) for _ in range(3)]
    temp_y = aimFunction(para,x)
    if temp_y - cur_best  < 0:
        Need=para
        cur_best=temp_y
        update_time+=1
        print("update!!",update_time,i)
result=[]
print(Need)
for i in x:
    result.append(myFuc(Need[0],Need[1],Need[2],i))
Draw([y,result],["real","predict"],x,"Real and predict",Colors=None)

备注

绘图函数，基于数学建模选拔赛代码整合封装得到

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.pyplot import MultipleLocator
plt.rcParams['font.sans-serif']=['SimHei'] #用来正常显示中文标签
plt.rcParams['axes.unicode_minus']=False #用来正常显示负号
#families=[ 'fantasy','Tahoma', 'monospace','Times New Roman']
def Draw(Lines,Names,index,Title,Colors=None):
    if Colors==None:
        Colors=['lightseagreen','orange','red','darkblue','green', ]
    plt.rcParams['figure.figsize'] = (12.0, 7.0)
    for i in range(len(Lines)):
        plt.plot(index,Lines[i],Colors[i],label=Names[i],linewidth=0.8)
    Title_font = {'family': 'fantasy','size': 20}
    plt.title(Title, Title_font)
    font1 = {'family': 'Tahoma','size': 14}
    plt.xlabel('X', font1)
    plt.ylabel('Y', font1,rotation=0)
    #plt.ylabel('Y', font1)
    plt.legend()
    x_len=len(index)/20
    x_major_locator = MultipleLocator(x_len)
    ax = plt.gca()
    ax.xaxis.set_major_locator(x_major_locator)
    plt.show()
#X = np.linspace(-np.pi, np.pi, 256, endpoint=True)
#C,S = np.cos(X), np.sin(X)
#Draw([C,S],['Sin','Cos'],X,"test")