2021-02-03美赛前MATLAB的学习笔记（全局优化、启发式算法、函数的写法）

全局优化

先给出一张表，并不表示我都会

算法	MATLAB求解器	作用
全局搜索	GlobalSearch	寻找全局最小
多起点搜索	MultiStart	寻找多个局部最小值
模式搜索	patternsearch	用模式搜索寻找函数最小值
遗传	Ga	用遗传算法寻找函数最小值
粒子群	particleswarm	用粒子群算法寻找函数的最小值
模拟退火	simulannealbnd	用模拟退火算法寻找函数最小值

这里再挖一个坑 NP问题
下面是一道经典的NP问题的题，TSP问题，也就是旅行推销员问题
精简的表述是：已给一个n个点的完全图，每条边都有一个长度，求总长度最短的经过每个顶点正好一次的封闭回路

对于这种问题，没有有效的算法，时间复杂度通常和暴力搜索等价
可以通过启发式算法得到近似最优解

遗传算法

遗传算法模仿大自然的遗传和自然选择，就像进化那样，后代种群比前代更适应环境，可以作为近似最优解。

遗传算法的参数

种群规模n，n不宜过大，会让收敛时间变长，n不宜过小，不能提供足够的采样点，一般在20~160
交叉概率$p_c$ 不宜过大，会把高适应的结构破坏，不宜太小，会让搜索停滞不前，一般在0.5~1
变异概率$p_m$ 太小，不会产生新的基因块，太大，会使遗传算法变成随机搜索，一般$p_m$取0.001~1
进化代数t 一般取100-1000，进化的选取，可以设定某种判定准则

染色体编码

常用的有二进制编码和浮点数编码，二进制编码直观但个体长度短时难以达到精度要求。由种群表现型组成的问题空间向种群基因型组成的编码空间的映射称作编码，反之是解码。

适应度函数

用来恒量个体优势，度量个体适应度的函数，适应度函数值越大的个体越好，一般来说，适应度函数是目标函数变换而成的。
假如,对于目标函数是最小化问题，可以建立如下适应度函数$F(x)$

$$F(x)=\{\begin{array}{c}{C}_{max}-f(x),f(x)<C_{max}\\ 0,f(x)>=C_{max}\end{array}$$
or
$$F(x)=\frac{1}{1+c+f(x)}$$

约束条件处理

罚函数法

不满足约束条件时，适应度函数减去或除以一个罚函数，难点在于度量是否满足约束条件的时候，降低了计算效率。

搜索空间限定法

剪枝，对搜索空间大小加以限制

遗传算子

遗传算子包括选择、交叉、变异，

选择

包括轮盘赌法，由适应值对应的概率决定被淘汰的概率
排序选择法
两两竞争法

交叉

单点、双点、多点、均匀交叉

变异

变异有基本位变异和均匀变异。基本位变异是对个体编码串中以变异概率、随机指定的某一位或某几位仅因座上的值做变异运算。均匀变异是分别用符合某一范围内均匀分布的随机数，以某一较小的概率来替换个体编码串中各个基因座上的原有基因值。

倒位

代码实现:非线性规划

现在要解决如下优化问题

$minf(x)=100(x_1^2-x_2{)}^2+(1-x_1{)}^2$
$$s.t.\{\begin{array}{c}{x}_1{x}_2+x_1-x_2+1.5<=0\\ 10-x_1{x}_2<=0\\ 0<=x_1<=1,0<=x_2<=13\end{array}$$

function y=simple_fitness(x)//目标函数
y=100*(x(1)^2-x(2))^2+(1-x(1))^2;

function [c,ceq]=simple_constraint(x)//约束函数
c=[1.5+x(1)*x(2)+x(1)-x(2)];10-x(1)*x(2)];
ceq=[];
%调用ga函数
ObjectiveFunction = @simple_fitness;
nvars=2;//变量数
LB=[0 0];//lower bound
UB=[1 13];//upper bound
ConstraintFunction = @simple_constraint;
[x，fval]=ga(ObjectiveFunction,nvars,[],[],[],[],LB,UB,ConstraintFunction)

让我们研究一下ga函数

x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCon,options)
or

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)

or

x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon)

省略以此类推

function

这函数的写法matlab就和其他语言好不同，直接在输出量里定义
这有必要单独列出来给你们瞧瞧

A,b,Aeq,beq

eq就是equal
A,b表示不等式 $a_1{x}_1+a_2{x}_2...<=b$
Aeq,beq表示等式

nonlcon

非线性约束，包含两个数组元素 ，c满足<=0,ceq满足=0

options

不会吧不会吧你想调参？
网址给你https://ww2.mathworks.cn/help/gads/ga.html?s_tid=srchtitle#mw_4a8bfdb9-7c4c-4302-8f47-d260b7a43e26

代码实现:TSP

加载可视化数据

load('usborder.mat','x','y','xx','yy');
plot(x,y,'Color','red');hold on;
cities=40;
locations=zeros(cities,2);
n=1;
while(n<=cities)
	xp=rand*1.5;
	yp=rand;
	if inpolygon(xp,yp,xx,yy)
		locations(n,1)=xp;
		locations(n,2)=yp;
		n=n+1;
	end
end
plot(locations(:,1),locations(:,2),'bo');
xlabel('城市横坐标x');
ylabel('城市纵坐标y');
hold on;

inpolygon & numel

确定位于多边形区域边缘内部或边缘上的点数。
多边形可以是点集

[in,on] = inpolygon(xq,yq,xv,yv);

in是在里面 on是在边缘上面

numel(xq(in))

numel统计数组个数

计算城市间距离

distances=zeros(cities);
for count1=1:cities
for count2=1:count1
	x1=locations(count1,1);
	y1=locations(count1,2);
	x2=locations(count2,1);
	y2=locations(count2,2);
	distances(count1,count2)=sqrt((x1-x2)^2+(y1-y2)^2);
	distances(count2,count1)=distances(count1,count2);
end
end	

定义目标函数和求解

FitnessFcn = @(x) traveling_salesman_fitness(x,distances);
  
my_plot = @(options,state,flag) traveling_salesman_plot(options, ...
    state,flag,locations);

options=optimoptions(@ga,'PopulationType','custom','InitialPopulationRange',[1;cities]);
options=optimoptions(options,
					'CreationFcn',@create_permutations,...
					'CrossoverFcn',@crossover_permutation,...
					'MutationFcn',@mutate_permutation,...
					'PlotFcn',my_plot,...
					'MaxStallGenerations',200,...
					'UseVectorized',true);
[x,fval,reason,output]=...
ga(FitnessFcn,cities,[],[],[],[],[],[],[],options);

其中traveling_salesman_fitness如下

function scores = traveling_salesman_fitness(x,distances)
%TRAVELING_SALESMAN_FITNESS  Custom fitness function for TSP. 
%   SCORES = TRAVELING_SALESMAN_FITNESS(X,DISTANCES) Calculate the fitness 
%   of an individual. The fitness is the total distance traveled for an
%   ordered set of cities in X. DISTANCE(A,B) is the distance from the city
%   A to the city B.

scores = zeros(size(x,1),1);
for j = 1:size(x,1)
    % here is where the special knowledge that the population is a cell
    % array is used. Normally, this would be pop(j,:);
    p = x{j}; 
    f = distances(p(end),p(1));
    for i = 2:length(p)
        f = f + distances(p(i-1),p(i));
    end
    scores(j) = f;
end

traveling_salesman_plot如下

function state = traveling_salesman_plot(options,state,flag,locations)
%   TRAVELING_SALESMAN_PLOT Custom plot function for traveling salesman.
%   STATE = TRAVELING_SALESMAN_PLOT(OPTIONS,STATE,FLAG,LOCATIONS) Plot city
%   LOCATIONS and connecting route between them. This function is specific
%   to the traveling salesman problem.

%   Copyright 2004-2006 The MathWorks, Inc.
persistent x y xx yy
if strcmpi(flag,'init')
  load('usborder.mat','x','y','xx','yy');
end
% plot(x,y,'Color','red');
% hold on;
[unused,i] = min(state.Score);
genotype = state.Population{i};

plot(locations(:,1),locations(:,2),'bo');
axis([-0.1 1.5 -0.2 1.2]);
hold on
plot(locations(genotype,1),locations(genotype,2));
xlabel('城市的横坐标x'); ylabel('城市的纵坐标y'); 
grid on
hold off

create_permutations如下

function pop = create_permutations(NVARS,FitnessFcn,options)
%CREATE_PERMUTATIONS Creates a population of permutations.
%   POP = CREATE_PERMUTATION(NVARS,FITNESSFCN,OPTIONS) creates a population
%  of permutations POP each with a length of NVARS. 
%
%   The arguments to the function are 
%     NVARS: Number of variables 
%     FITNESSFCN: Fitness function 
%     OPTIONS: Options structure used by the GA

%   Copyright 2004-2007 The MathWorks, Inc.

totalPopulationSize = sum(options.PopulationSize);
n = NVARS;
pop = cell(totalPopulationSize,1);
for i = 1:totalPopulationSize
    pop{i} = randperm(n); 
end

crossover_permutation如下

function xoverKids  = crossover_permutation(parents,options,NVARS, ...
    FitnessFcn,thisScore,thisPopulation)
%   CROSSOVER_PERMUTATION Custom crossover function for traveling salesman.
%   XOVERKIDS = CROSSOVER_PERMUTATION(PARENTS,OPTIONS,NVARS, ...
%   FITNESSFCN,THISSCORE,THISPOPULATION) crossovers PARENTS to produce
%   the children XOVERKIDS.
%
%   The arguments to the function are 
%     PARENTS: Parents chosen by the selection function
%     OPTIONS: Options created from OPTIMOPTIONS
%     NVARS: Number of variables 
%     FITNESSFCN: Fitness function 
%     STATE: State structure used by the GA solver 
%     THISSCORE: Vector of scores of the current population 
%     THISPOPULATION: Matrix of individuals in the current population

%   Copyright 2004-2015 The MathWorks, Inc. 

nKids = length(parents)/2;
xoverKids = cell(nKids,1); % Normally zeros(nKids,NVARS);
index = 1;

for i=1:nKids
    % here is where the special knowledge that the population is a cell
    % array is used. Normally, this would be thisPopulation(parents(index),:);
    parent = thisPopulation{parents(index)};
    index = index + 2;

    % Flip a section of parent1.
    p1 = ceil((length(parent) -1) * rand);
    p2 = p1 + ceil((length(parent) - p1- 1) * rand);
    child = parent;
    child(p1:p2) = fliplr(child(p1:p2));
    xoverKids{i} = child; % Normally, xoverKids(i,:);
end

mutate_permutation 如下

function mutationChildren = mutate_permutation(parents ,options,NVARS, ...
    FitnessFcn, state, thisScore,thisPopulation,mutationRate)
%   MUTATE_PERMUTATION Custom mutation function for traveling salesman.
%   MUTATIONCHILDREN = MUTATE_PERMUTATION(PARENTS,OPTIONS,NVARS, ...
%   FITNESSFCN,STATE,THISSCORE,THISPOPULATION,MUTATIONRATE) mutate the
%   PARENTS to produce mutated children MUTATIONCHILDREN.
%
%   The arguments to the function are 
%     PARENTS: Parents chosen by the selection function
%     OPTIONS: Options created from OPTIMOPTIONS
%     NVARS: Number of variables 
%     FITNESSFCN: Fitness function 
%     STATE: State structure used by the GA solver 
%     THISSCORE: Vector of scores of the current population 
%     THISPOPULATION: Matrix of individuals in the current population
%     MUTATIONRATE: Rate of mutation

%   Copyright 2004-2015 The MathWorks, Inc.

% Here we swap two elements of the permutation
mutationChildren = cell(length(parents),1);% Normally zeros(length(parents),NVARS);
for i=1:length(parents)
    parent = thisPopulation{parents(i)}; % Normally thisPopulation(parents(i),:)
    p = ceil(length(parent) * rand(1,2));
    child = parent;
    child(p(1)) = parent(p(2));
    child(p(2)) = parent(p(1));
    mutationChildren{i} = child; % Normally mutationChildren(i,:)
end

再研究ga

默认的目标函数输入是nvar个参数，如果option里的‘UseVectorized’是True的话，支持输入nvar*pop的向量，pop是会变的，所以要求输入能变

模拟退火

（挖坑挖坑Orz明天真有工作了，完了）