matlab粒子群算法

求方程组

vpasolve函数

用于求解符号方程组的函数

clear; clc
syms x1 x2 x3
eqn =  [abs(x1+x2)-abs(x3)  == 0, x1*x2*x3+18 == 0, x1^2*x2+3*x3 == 0]
[x1, x2, x3] = vpasolve(eqn, [x1, x2, x3], [0 0 0])  
[x1, x2, x3] = vpasolve(eqn, [x1, x2, x3], [1 1 1])  % 换一个初始值试试

注释：

1.定义符号变量x1 x2 x3，

syms x1 x2 x3

2.注意符号方程要写；两个等号：==

abs(x1+x2)-abs(x3)  == 0

3.设置初始值==[0 0 0]==

[x1, x2, x3] = vpasolve(eqn, [x1, x2, x3], [0 0 0])  

fsolve函数

x0 = [0,0,0];  % 初始值
x = fsolve(@my_fun,x0)
% 换一个初始值试试
x0 = [1,1,1];  % 初始值
format long g  % 显示更多的小数位数
x = fsolve(@my_fun,x0)

注释：

把约束条件写到==@my_fun==文件中，

粒子群算法求解

需要先构造目标函数：

函数代码：

function f = Obj_fun(x)
    f1= abs(x(1)+x(2))-abs(x(3)) ;
    f2 = x(1) * x(2) * x(3) + 18;
    f3= x(1)^2 * x(2) + 3*x(3);
    f = abs(f1) + abs(f2) + abs(f3);
end

注释：
1.

代码：

clear; clc
narvs = 3;
% 使用粒子群算法，不需要指定初始值，只需要给定一个搜索的范围
lb = -10*ones(1,3);  ub = 10*ones(1,3);  
options = optimoptions('particleswarm','FunctionTolerance',1e-12,'MaxStallIterations',100,'MaxIterations',20000,'SwarmSize',100);
[x, fval] = particleswarm(@Obj_fun,narvs,lb,ub,options) 

注释：
1.narvs = 3;设置变量个数。
2.lb = -10*ones(1,3); ub = 10*ones(1,3); 设置变量最大值和最小值。

3.可能有多个解：多次求解尝试。

多元函数拟合：

f竖线后的参数即为未知的待估参数。

粒子群算法求解多元函数拟合

求解微分方差代码：

clear;clc
load mydata.mat  % 导入数据（共三列，分别是S,I,R的数量）
n = size(mydata,1);  % 一共有多少行数据
true_s = mydata(:,1);   
true_i = mydata(:,2);
true_r = mydata(:,3);
figure(1) 
plot(1:n,true_s,'r-',1:n,true_i,'b-',1:n,true_r,'g-')  % S的数量太大了
legend('S','I','R')
plot(1:n,true_i,'b-',1:n,true_r,'g-') % 单独画出真实的I和R的数量
hold on  % 等会接着在这个图形上面画图

% 随便先给一组参数来计算微分方程的数值解
beta = 0.1;   % 易感染者与已感染者接触且被传染的强度
gamma = 0.01;  % 康复率
[t,p]=ode45(@(t,x) SIR_fun(t,x,beta,gamma), [1:n],[true_s(1) true_i(1) true_r(1)]);
% p就是计算的数值解，它有三列，分别是S I R的数量
p = round(p);  % 对p进行四舍五入(人数为整数)
% 注意，四舍五入后有可能出现总人数之和不为10000的情况，例如为9999或10001，但是这种误差可以忽略不计。
plot(1:n,p(:,2),'b--',1:n,p(:,3),'g--') 
legend('I','R','拟合的I','拟合的R')
hold off  % 和hold on对应
% 计算残差平方和
sse = sum((true_s - p(:,1)).^2  + (true_i -  p(:,2)).^2  + (true_r - p(:,3)).^2)

解释：

[t,p]=ode45(@(t,x) SIR_fun(t,x,beta,gamma), [1:n],[true_s(1) true_i(1) true_r(1)]);

ode45函数：（求微分方程数值解）

ode45函数
第一个参数：$dx=f(t)$
第二个参数：是自变量的取值范围。
第三个参数：微分方程因变量初始值。
函数代码：

function dx=SIR_fun(t,x,beta,gamma)   
    % beta: 易感染者与已感染者接触且被传染的强度
    % gamma: 康复率
    dx = zeros(3,1);  % x(1)表示S  x(2)表示I  x(3)表示R
    C = x(1)+x(2);  % 传染病系统中的有效人群,也就是课件中的N' = S+I
    dx(1) = - beta*x(1)*x(2)/C;  
    dx(2) = beta*x(1)*x(2)/C - gamma*x(2);
    dx(3) = gamma*x(2); 
end

网格搜索BETA与GAMMA

粒子群算法工具箱函数

options = optimoptions('particleswarm','Display','iter','SwarmSize',100,'PlotFcn','pswplotbestf');
[h, fval] = particleswarm(@Obj_fun,2,lb,ub,options)  % h就是要优化的参数，fval是目标函数值

注释：
1.@Obj_fun为目标函数。
2.2求解参数个数。
3.lb,ub求解参数的上下界。
4.options选项，

options = optimoptions('particleswarm','Display','iter','SwarmSize',100,'PlotFcn','pswplotbestf');

Display为在命令行窗口展示迭代的过程。
'SwarmSize',100设置粒子群迭代的大小，
'PlotFcn','pswplotbestf'绘制最优的适应度变化图。

生成结果如图：
我们重点关注第三列。（最佳的适应度变化图）（即最小的目标函数）
本题为SSE（残差平方和）

function dx=SIR_fun(t,x,beta,gamma)   
     % beta: 易感染者与已感染者接触且被传染的强度
     % gamma: 康复率
    dx = zeros(3,1);  % x(1)表示S  x(2)表示I  x(3)表示R
    C = x(1)+x(2);  % 传染病系统中的有效人群,也就是课件中的N' = S+I
    dx(1) = - beta*x(1)*x(2)/C;  
    dx(2) = beta*x(1)*x(2)/C - gamma*x(2);
    dx(3) = gamma*x(2); 
end

x为一个列向量，分别代表S，I，R三个因变量

目标函数：

function f = Obj_fun(h)  
    % 注意，h就是我们要估计的目标函数里面的参数
    % h的长度就是我们待估参数的个数
    global true_s  true_i  true_r n   % 在子函数中使用全局变量前也需要声明下
    beta = h(1);  % 易感染者与已感染者接触且被传染的强度
    gamma = h(2);  % 康复率
    [t,p]=ode45(@(t,x) SIR_fun(t,x,beta,gamma), [1:n],[true_s(1) true_i(1) true_r(1)]);   
    p = round(p);  % 对p进行四舍五入(人数为整数)
    % p的第一列是易感染者S的数量，p的第二列是患者I的数量，p的第三列是康复者R的数量
    f = sum((true_s - p(:,1)).^2  + (true_i -  p(:,2)).^2  + (true_r - p(:,3)).^2);
end

beta和gamma是随时间变换的函数的拟合

有四个变量需要被拟合，
目标函数长度应该为4，