目前我的这个代码第八行报错,基于pso的电池离线辨识的一个代码,请大家看一下
%% 程序初始化 √
% x = [0.1 0.155 0.2 15 15]
clc;
clear;
close all;
%% 初始参数 √
data=xlsread('D:\仿真真值数据.xlsx');
data=data( 1:43932 ,:);%1-0.9
I=data(:,2);%实际电流
Ut=data(:,3);%实际电压
SOC=data(:,4);%真实SOC值----安时积分法
OCVp=[1.736 -1.306 -10.7 22.41 -16.73 5.96 2.738];%OCV_SOC曲线拟合系数
OCV=polyval(OCVp,SOC);%开路电压值
%% 优化相关参数设置 √
N = 100; % 初始种群个数100
d = 5; % 空间维数
ger = 200; % 最大迭代次数200
Rlimit = [0.1, 0.2]; % 设置x参数限制
Climit = [12, 17]; % 设置x参数限制
vRlimit = [-0.1,0.1]; % 设置速度限制
vClimit = [-1,1]; % 设置速度限制
% w = 0.8; % 惯性权重
wmax = 0.9; % 惯性权重
wmin = 0.4; % 惯性权重
c1 = 1.5; % 自我学习因子
c2 = 1.49; % 群体学习因子
%% 初始化种群 √
pop = zeros(N,d);
for i = 1:N % 对每个个体进行位置的初始化,需要初始化的参数为x和y
pop(i,1) = Rlimit(1) + (Rlimit(2) - Rlimit(1)) * rand; % 初始化参数x
pop(i,2) = Rlimit(1) + (Rlimit(2) - Rlimit(1)) * rand; % 初始化参数x
pop(i,3) = Rlimit(1) + (Rlimit(2) - Rlimit(1)) * rand; % 初始化参数x
pop(i,4) = Climit(1) + (Climit(2) - Climit(1)) * rand; % 初始化参数x
pop(i,5) = Climit(1) + (Climit(2) - Climit(1)) * rand; % 初始化参数x
end
v = rand(N, d); % 初始种群的速度,直接在0~1内随机生成一个N*d的数据矩阵
local_best = pop; % 每个个体的历史最佳位置,N*d矩阵
global_best = zeros(1, d); % 种群的历史最佳位置,初始化种群的最佳位置为原点,1*d矩阵
f_local_best = zeros(N, 1); % 每个个体的历史最佳适应度,N*1矩阵
f_global_best = 0; % 种群历史最佳适应度
%% 更新种群 √
iter = 1; % 当前迭代次数
record = zeros(ger, 1); % 记录器,用于观察收敛程度
while iter <= ger % 当当前迭代次数小于最大迭代次数时
w=wmax-((wmax-wmin)*(iter^2))/(ger^2);
parfor j = 1:N
f_pop(j) = fitness(pop(j,:),I,Ut,SOC,OCV) ;% f_pop指代个体当前适应度,每个个体的x和y存储在pop矩阵中
end
iter;
for i = 1:N % 更新每个个体的历史最佳适应度和历史最佳位置
if f_local_best(i) < f_pop(i)
f_local_best(i) = f_pop(i); % 当个体适应度大于个体最佳适应度时,更新个体历史最佳适应度
local_best(i,:) = pop(i,:); % 更新第i个个体的历史最佳位置
end
end
if f_global_best < max(f_local_best)
[f_global_best, nmax] = max(f_local_best);
global_best = pop(nmax, :); % 更新群体历史最佳位置
end
v = v * w + c1 * rand * (local_best - pop) + c2 * rand * (repmat(global_best, N, 1) - pop); % 速度更新公式,!!注意使用repmat命令后矩阵的维数!!
vR_up = find(v(:,1:3) > vRlimit(2)); % v_up指代速度矩阵v中,即N*d维矩阵中数值大于v上限的元素的位置
vR_low = find(v(:,1:3) < vRlimit(1)); % v_low指代速度矩阵v中,即N*d维矩阵中数值小于v下限的元素的位置
v(vR_up) = vRlimit(2); % 将速度矩阵v中大于v上限的元素数值变为v上限
v(vR_low) = vRlimit(1); % 将速度矩阵v中小于v下限的元素数值变为v下限
vC_up = find(v(:,4:5) > vClimit(2)); % v_up指代速度矩阵v中,即N*d维矩阵中数值大于v上限的元素的位置
vC_low = find(v(:,4:5) < vClimit(1)); % v_low指代速度矩阵v中,即N*d维矩阵中数值小于v下限的元素的位置
v(vC_up) = vClimit(2); % 将速度矩阵v中大于v上限的元素数值变为v上限
v(vC_low) = vClimit(1); % 将速度矩阵v中小于v下限的元素数值变为v下限
pop = pop + v;
pop_R_up = find(pop(:,1:3) > Rlimit(2)); % pop_x_up指代位置矩阵pop第一列,即每个个体x取值中大于x上限的位置
pop_R_low = find(pop(:,1:3) < Rlimit(1)); % pop_x_low指代位置矩阵pop第一列,即每个个体x取值中小于x下限的位置
pop(pop_R_up) = Rlimit(2); % 将位置矩阵pop第一列中大于x上限的取值变为x上限
pop(pop_R_low) = Rlimit(1); % 将位置矩阵pop第一列中小于x下限的取值变为x下限
pop_C_up = find(pop(:,4:5) > Climit(2)); % pop_x_up指代位置矩阵pop第一列,即每个个体x取值中大于x上限的位置
pop_C_low = find(pop(:,4:5) < Climit(1)); % pop_x_low指代位置矩阵pop第一列,即每个个体x取值中小于x下限的位置
pop(pop_C_up) = Climit(2); % 将位置矩阵pop第一列中大于x上限的取值变为x上限
pop(pop_C_low) = Climit(1); % 将位置矩阵pop第一列中小于x下限的取值变为x下限
record(iter) = f_global_best; %最大值记录
iter = iter+1;
end % 一次迭代完成
figure(1)
plot(record);
title('收敛过程')
xlabel('迭代次数')
ylabel('适应度')
hold on
disp(['最大值:',num2str(f_global_best)]);
disp(['R0 = ',num2str(global_best(1,1))]);
disp(['R1 = ',num2str(global_best(1,2))]);
disp(['R2 = ',num2str(global_best(1,3))]);
disp(['C1 = ',num2str(global_best(1,4))]);
disp(['C2 = ',num2str(global_best(1,5))]);
%% 求最小端电压误差
R0=global_best(1,1);
R1=global_best(1,2);
R2=global_best(1,3);
C1=global_best(1,4);
C2=global_best(1,5);
T=size(I,1);%电流的行数
Ut_predict=zeros(T,1);%预测电压
Ut_error=zeros(T,1);%预测电压误差
T_s=1;%采样时间
U1=0;%RC电压
U2=0;%RC电压
Ut_error_mse=0;
for k=1:T
U1=exp(-T_s/(R1*C1))*U1+R1*(1-exp(-T_s/(R1*C1)))*I(k);%极化电压
U2=exp(-T_s/(R2*C2))*U2+R2*(1-exp(-T_s/(R2*C2)))*I(k);%极化电压
Ut_predict(k)=OCV(k)+U1+U2+R0*I(k);%端电压预测值
Ut_error(k)=Ut(k)-Ut_predict(k);%端电压预测误差
Ut_error_mse=Ut_error_mse+(Ut_error(k))^2;%0-k时刻误差平方和(没有乘以时间的)
end
Ut_error;
%% 画图
%端电压预测图
figure(2)
hold on
plot(Ut,'k')
plot(Ut_predict,'r--')
legend('真实电压','预测电压')
xlabel('时间(s)')
ylabel('端电压(V)')
title('端电压图')
%端电压误差图
figure(3)
plot(Ut_error)
xlabel('时间(s)')
ylabel('误差(V)')
title('端电压误差图')