基于纳什谈判理论的风–光–氢多主体能源系统合作运行方法（含matlab代码）

maxIter=50;%最大迭代次数
tolerant=1e-5;%收敛精度
iter=1;%迭代次数
Ben_Store=[];%历史目标函数
toler1=[];%残差1，风电主体
toler2=[];%残差2，光伏主体
P_pv2h_2=zeros(maxIter+1,24);P_wt2h_2=zeros(maxIter+1,24);
value_P_wt2h_2=zeros(maxIter+1,24);value_P_pv2h_2=zeros(maxIter+1,24);
%% 迭代
while 1
    if iter==maxIter  %限制迭代次数
       disp('迭代不收敛,参数有误');
       break; 
    end 
    display(['迭代还未收敛,当前迭代第 ', num2str(iter),' 次']);
    if iter==1  %第一次求解比较特殊，要给一致性变量P_wt2h_1、P_pv2h_1赋初值
        [P_pv2h_2(2,:),P_wt2h_2(2,:),Obj_dzq(iter)]=Fun_dzq(value_P_pv2h_2(iter,:),value_P_wt2h_2(iter,:),rho_wt,rho_pv,lambda_wt,lambda_pv);
        [value_P_pv2h_2(2,:),Obj_gf(iter)]=Fun_gf(P_pv2h_2(2,:),rho_pv,lambda_pv);
        [value_P_wt2h_2(2,:),Obj_fd(iter)]=Fun_fd(P_wt2h_2(2,:),rho_wt,lambda_wt);
        lambda_wt=lambda_wt+rho_wt*(P_wt2h_2(2,:)-value_P_wt2h_2(2,:));
        lambda_pv=lambda_pv+rho_pv*(P_pv2h_2(2,:)-value_P_pv2h_2(2,:));
    else
        [P_pv2h_2(iter+1,:),P_wt2h_2(iter+1,:),Obj_dzq(iter)]=Fun_dzq(value_P_pv2h_2(iter,:),value_P_wt2h_2(iter,:),rho_wt,rho_pv,lambda_wt,lambda_pv);
        [value_P_pv2h_2(iter+1,:),Obj_gf(iter)]=Fun_gf(P_pv2h_2(iter+1,:),rho_pv,lambda_pv);
        [value_P_wt2h_2(iter+1,:),Obj_fd(iter)]=Fun_fd(P_wt2h_2(iter+1,:),rho_wt,lambda_wt);
        lambda_wt=lambda_wt+rho_wt*(P_wt2h_2(iter+1,:)-value_P_wt2h_2(iter+1,:));
        lambda_pv=lambda_pv+rho_pv*(P_pv2h_2(iter+1,:)-value_P_pv2h_2(iter+1,:));
    end
    %保存历史数据
    Ben_Store=[Ben_Store,[Obj_dzq(iter);Obj_gf(iter);Obj_fd(iter)]];
    %残差计算,公式（36）应该修改一下,是向量的二范数再平方
    toler1=[toler1,norm(P_wt2h_2(iter+1,:)-value_P_wt2h_2(iter+1,:))^2];%保存残差1 
    toler2=[toler2,norm(P_pv2h_2(iter+1,:)-value_P_pv2h_2(iter+1,:))^2];%保存残差2                
    %判断收敛条件
    if toler1(iter)<=tolerant && toler2(iter)<=tolerant
       display(['迭代收敛,在第 ', num2str(iter),' 次收敛']);
       break; 
    end
    iter=iter+1;
end
%% 画图
figure(1)
plot(-Ben_Store(1,:),'b-o','LineWidth',1.5);
xlabel('迭代次数');
ylabel('成本/元');
title('电制氢主体分布式迭代情况');
;
figure(2)
plot(-Ben_Store(2,:),'r-o','LineWidth',1.5);
xlabel('迭代次数');
ylabel('成本/元');
title('光伏主体分布式迭代情况');
;
figure(3)
plot(-Ben_Store(3,:),'k-o','LineWidth',1.5);
xlabel('迭代次数');
ylabel('成本/元');
title('风力主体分布式迭代情况');
;
%
Party=-Ben_Store(1,:)-Ben_Store(2,:)-Ben_Store(3,:);
figure(4)
plot(Party,'k-o','LineWidth',1.5);
xlabel('迭代次数');
ylabel('成本/元');
title('风光联盟总效益值');