118号资源-源程序：论文可在知网下载《基于纳什谈判理论的风–光–氢多主体能源系统 合作运行方法》本人博客有解读

clc clear all close all %% 0.部分变量名称对照表 % w_H 问题1求得的 电制氢主体最优解 value=-4909.15746156799 % w_WT 问题1求得的 风电主体最优解 value=6709.93831829112 % w_PV 问题1求得的 光伏主体最优解 value=8428.84347521765 % P_wt2h 问题1求得的 风-氢之间的最优期望交易电量 % P_pv2h 问题1求得的 风-光之间的最优期望交易电量 % U_0_H2 谈判破裂点 value=-15798 % U_0_wt 谈判破裂点 value=12598 % U_0_pv 谈判破裂点 value=8579 %% 1.常量赋值 P_wt2h=[0,0,0,0,0,0,0.000118965046880248,1047.26731848920,1047.26724108527,1047.26712789578,1047.26702135156,1047.26691403564,1047.26680723930,1047.26670219614,824.427467809682,786.259531081652,1047.26638935748,1047.26621375800,1047.26609919524,1047.26624366009,1047.26610449827,1047.26598842897,0,0]; P_pv2h=[0,0,0,0,0,0,0,47.2692475009289,47.2691745816221,47.2691054587244,47.2690389348969,47.2689731137088,47.2689067230860,47.2688424477525,47.2701311975193,47.2703042915214,47.2686462019993,47.2685518537297,47.2684560245753,0,0,0,0,0]; %% 2.变量声明 % P_wt_1=sdpvar(1,24);% 光伏主体： 风-氢交易电价 % P_pv_1=sdpvar(1,24);% 风电主体： 风-光交易电价 % P_wt_2=sdpvar(1,24);% 电制氢主体：风-氢交易电价 % P_pv_2=sdpvar(1,24);% 电制氢主体：风-光交易电价 %% 3.ADMM迭代参数设置 rho_wt=1;%惩罚因子 rho_pv=1;%惩罚因子 lambda_wt=0*ones(1,24);%风主体拉格朗日乘子 lambda_pv=0*ones(1,24);%光主体拉格朗日乘子 maxIter=50;%最大迭代次数 tolerant=1e-4;%收敛精度 iter=1;%迭代次数 Ben_Store=[];%历史目标函数 toler1=[];%残差1，风电主体 toler2=[];%残差2，光伏主体 P_pv_2=zeros(maxIter+1,24);P_wt_2=zeros(maxIter+1,24); value_P_wt_2=zeros(maxIter+1,24);value_P_pv_2=zeros(maxIter+1,24); %% 4.迭代 while 1 if iter==maxIter %限制迭代次数 disp('迭代不收敛,参数有误'); break; end display(['迭代还未收敛,当前迭代第 ', num2str(iter),' 次']); if iter==1 %第一次求解比较特殊，要给一致性变量P_wt_1、P_pv_1赋初值 [P_pv_2(2,:),P_wt_2(2,:),Obj_dzq(iter)]=Fun2_dzq(value_P_wt_2(iter,:),value_P_pv_2(iter,:),rho_wt,rho_pv,lambda_wt,lambda_pv); [value_P_pv_2(2,:),Obj_gf(iter)]=Fun2_gf(P_pv_2(2,:),rho_pv,lambda_pv); [value_P_wt_2(2,:),Obj_fd(iter)]=Fun2_fd(P_wt_2(2,:),rho_wt,lambda_wt); lambda_wt=lambda_wt+rho_wt*(P_wt_2(2,:)-value_P_wt_2(2,:)); lambda_pv=lambda_pv+rho_pv*(P_pv_2(2,:)-value_P_pv_2(2,:)); else [P_pv_2(iter+1,:),P_wt_2(iter+1,:),Obj_dzq(iter)]=Fun2_dzq(value_P_wt_2(iter,:),value_P_pv_2(iter,:),rho_wt,rho_pv,lambda_wt,lambda_pv); [value_P_pv_2(iter+1,:),Obj_gf(iter)]=Fun2_gf(P_pv_2(iter+1,:),rho_pv,lambda_pv); [value_P_wt_2(iter+1,:),Obj_fd(iter)]=Fun2_fd(P_wt_2(iter+1,:),rho_wt,lambda_wt); lambda_wt=lambda_wt+rho_wt*(P_wt_2(iter+1,:)-value_P_wt_2(iter+1,:)); lambda_pv=lambda_pv+rho_pv*(P_pv_2(iter+1,:)-value_P_pv_2(iter+1,:)); end %%%%%%%%%%%%保存历史数据 Ben_Store=[Ben_Store,[Obj_dzq(iter);,Obj_gf(iter);,Obj_fd(iter)]]; %残差计算，公式（36）应该修改一下，是 向量的二范数再平方 toler1=[toler1,norm(P_wt_2(iter+1,:)-value_P_wt_2(iter+1,:))^2];%保存残差1 toler2=[toler2,norm(P_pv_2(iter+1,:)-value_P_pv_2(iter+1,:))^2];%保存残差2 %%%%%%%%%%%%判断收敛条件 if toler1(iter)<=tolerant && toler2(iter)<=tolerant display(['迭代收敛,在第 ', num2str(iter),' 次收敛']); break; end iter=iter+1; end %% 画图 figure plot(toler1,'b-o','LineWidth',1.5); xlabel('迭代次数'); ylabel(' '); title('残差1'); box off figure plot(toler2,'r-o','LineWidth',1.5); xlabel('迭代次数'); ylabel(' '); title('残差2'); box off figure plot(P_wt_2(iter,:),'k-o','LineWidth',1.5); xlabel('迭代次数'); ylabel('电价/（元/（kW·h））'); title('风-氢交易电价'); box off figure plot(P_pv_2(iter,:),'k-o','LineWidth',1.5); xlabel('迭代次数'); ylabel('电价/（元/（kW·h））'); title('光-氢交易电价'); box off % ops=sdpsettings('solver','IPOPT','verbose',0,'usex0',0); % ops.cplex.mip.tolerances.mipgap=1e-6; % result=solvesdp(C,Obj,ops); % mosek_path = 'C:\Program Files\mosek\9.2\toolbox\R2015a'; % addpath(mosek_path); % [x,fval,exitflag,output,lambda] = quadprog(C,Obj); % figure(1) % plot(double(P_wt_1)) % figure(2) % plot(double(P_pv_1))