clear
clc
close all;
load('srdpshuju')
load('dfmat')
T0=10000; %初始温度
T=T0;
r=0.975; %衰减指数(指数式衰减)
N=100; %解的规模(并行)
M=20; %Markov链长度
b=10; %概率函数中的系数
%K=1.3806488*10^-23; %玻尔兹曼常数
a0=8.8542*10^-12; %真空中介电常数
u0=4*pi*10^-7; %真空中磁导率
c=3*10^8; %光速
F=3*10^8; %主频
fl=10; %频率下限
fh=40; %频率上限
z=zeros(1,N);
s=struct('a1',z,'a2',z,'a3',z,'b1',z,'b2',z,'t0',z,'T',z); % a是介电常数,b是电导率
ss=struct('a1',z,'a2',z,'a3',z,'b1',z,'b2',z,'t0',z,'T',z); %搜索解
Z=0;
sz=struct('a1',Z,'a2',Z,'a3',Z,'b1',Z,'b2',Z,'t0',Z,'T',Z); %找100个中最优
sy=struct('a1',Z,'a2',Z,'a3',Z,'b1',Z,'b2',Z,'t0',Z,'T',Z); %找每次中最优
%下面采用rand函数对每个基因赋值
rng shuffle %时变种子,也可用rand('state',sum(clock))
s.a1=5+10*rand(1,N); %a,b根据实际情况要调整
s.a2=20+10*rand(1,N);
s.a3=5+10*rand(1,N);
s.b1=0.01+0.02*rand(1,N);
s.b2=0.005+0.01*rand(1,N);
s.t0=6.0*10^-9+4*10^-9*rand(1,N);
s.T=0.5*10^-9+1*10^-9*rand(1,N);
% rng shuffle %时变种子,也可用rand('state',sum(clock))
% s.a1=1+14*rand(1,N); %a,b根据实际情况要调整
% s.a2=15+18*rand(1,N);
% s.a3=3+10*rand(1,N);
% s.b1=0.0+0.04*rand(1,N);
% s.b2=0.0+0.02*rand(1,N);
% s.t0=5.0*10^-9+6*10^-9*rand(1,N);
% s.T=0*10^-9+5*10^-9*rand(1,N);
Q1=zeros(1,N);
for k=1:N
Q1(k)=0;
for f=fl:fh
v1=c/sqrt(s.a1(k));
v2=c/sqrt(s.a2(k));
z0=0; %水平层厚度
z1=v1*s.t0(k)/2;
z2=v2*s.T(k)/2;
A0=0; %衰减常数
m1=(u0*s.a1(k)*a0)/2;
n1=sqrt(1+(s.b1(k)/(2*pi*(f-1)*df*s.a1(k)*a0))^2)-1;
A1=2*pi*(f-1)*df*sqrt(m1*n1); %
m2=(u0*s.a2(k)*a0)/2;
n2=sqrt(1+(s.b2(k)/(2*pi*(f-1)*df*s.a2(k)*a0))^2)-1;
A2=2*pi*(f-1)*df*sqrt(m2*n2); %
r01=(sqrt(a0)-sqrt(s.a1(k)*a0))/(sqrt(a0)+sqrt(s.a1(k)*a0)); %反射系数
r12=(sqrt(s.a1(k)*a0)-sqrt(s.a2(k)*a0))/(sqrt(s.a1(k)*a0)+sqrt(s.a2(k)*a0));
r23=(sqrt(s.a2(k)*a0)-sqrt(s.a3(k)*a0))/(sqrt(s.a2(k)*a0)+sqrt(s.a3(k)*a0));
r1=r12*(1-r01^2)*exp(-2*A1*z1)*exp(-A0*z0); %广义反射系数
r2=r23*(1-r01^2)*(1-r12^2)*exp(-2*A1*z1)*exp(-2*A2*z2)*exp(-A0*z0);
re=(r1+r2)/2;ro=(r1-r2)/2; %偶奇分量
rf=ref_fft(f); % 反射系数频域序列
q1=2*re*cos(pi*(f-1)*df*s.T(k))*cos(pi*(f-1)*df*(2*s.t0(k)+s.T(k)));
q2=2*ro*sin(pi*(f-1)*df*s.T(k))*sin(pi*(f-1)*df*(2*s.t0(k)+s.T(k)));
q3=2*re*cos(pi*(f-1)*df*s.T(k))*sin(pi*(f-1)*df*(2*s.t0(k)+s.T(k)));
q4=2*ro*sin(pi*(f-1)*df*s.T(k))*cos(pi*(f-1)*df*(2*s.t0(k)+s.T(k)));
q=abs(real(rf)-q1-q2)+abs(imag(rf)+q3-q4);
Q1(k)=Q1(k)+q; %目标函数值(适应度函数)
% v1=c/sqrt(s.a1(k));
% v2=c/sqrt(s.a2(k));
% z0=0; %水平层厚度
% z1=v1*s.t0(k)/2;
% z2=v2*s.T(k)/2;
%
% A0=0; %衰减常数
% m1=(u0*s.a1(k)*a0)/2;
% n1=sqrt(1+(s.b1(k)/(2*pi*(f)*df*s.a1(k)*a0))^2)-1;
% A1=2*pi*(f)*df*sqrt(m1*n1); %
% m2=(u0*s.a2(k)*a0)/2;
% n2=sqrt(1+(s.b2(k)/(2*pi*(f)*df*s.a2(k)*a0))^2)-1;
% A2=2*pi*(f)*df*sqrt(m2*n2); %
%
% r01=(sqrt(a0)-sqrt(s.a1(k)*a0))/(sqrt(a0)+sqrt(s.a1(k)*a0)); %反射系数
% r12=(sqrt(s.a1(k)*a0)-sqrt(s.a2(k)*a0))/(sqrt(s.a1(k)*a0)+sqrt(s.a2(k)*a0));
% r23=(sqrt(s.a2(k)*a0)-sqrt(s.a3(k)*a0))/(sqrt(s.a2(k)*a0)+sqrt(s.a3(k)*a0));
%
% r1=r12*(1-r01^2)*exp(-2*A1*z1)*exp(-A0*z0); %广义反射系数
% r2=r23*(1-r01^2)*(1-r12^2)*exp(-2*A1*z1)*exp(-2*A2*z2)*exp(-A0*z0);
% re=(r1+r2)/2;ro=(r1-r2)/2; %偶奇分量
%
% rf=ref_fft(f); % 反射系数频域序列
% q1=2*re*cos(pi*(f)*df*s.T(k))*cos(pi*(f)*df*(2*s.t0(k)+s.T(k)));
% q2=2*ro*sin(pi*(f)*df*s.T(k))*sin(pi*(f)*df*(2*s.t0(k)+s.T(k)));
% q3=2*re*cos(pi*(f)*df*s.T(k))*sin(pi*(f)*df*(2*s.t0(k)+s.T(k)));
% q4=2*ro*sin(pi*(f)*df*s.T(k))*cos(pi*(f)*df*(2*s.t0(k)+s.T(k)));
% q=abs(real(rf)-q1-q2)+abs(imag(rf)+q3-q4);
% Q1(k)=Q1(k)+q; %目标函数值(适应度函数)
end
end
d1=0.1;d2=0.2;d3=0.1;d4=0.001;d5=0.0005;d6=0.1*10^-9;d7=0.01*10^-9;
sign=randi([-1,1]);
ss.a1=s.a1+d1*sign*rand(1,N);
sign=randi([-1,1]);
ss.a2=s.a2+d2*sign*rand(1,N);
sign=randi([-1,1]);
ss.a3=s.a3+d3*sign*rand(1,N);
sign=randi([-1,1]);
ss.b1=s.b1+d4*sign*rand(1,N);
sign=randi([-1,1]);
ss.b2=s.b2+d5*sign*rand(1,N);
sign=randi([-1,1]);
ss.t0=s.t0+d6*sign*rand(1,N);
sign=randi([-1,1]);
ss.T=s.T+d7*sign*rand(1,N);
Q2=zeros(1,N);
for k=1:N
Q2(k)=0;
for f=fl:fh
v1=c/sqrt(ss.a1(k));
v2=c/sqrt(ss.a2(k));
z0=0; %水平层厚度
z1=v1*ss.t0(k)/2;
z2=v2*ss.T(k)/2;
A0=0; %衰减常数
m1=(u0*ss.a1(k)*a0)/2;
n1=sqrt(1+(ss.b1(k)/(2*pi*(f-1)*df*ss.a1(k)*a0))^2)-1;
A1=2*pi*(f-1)*df*sqrt(m1*n1); %
m2=(u0*ss.a2(k)*a0)/2;
n2=sqrt(1+(ss.b2(k)/(2*pi*(f-1)*df*ss.a2(k)*a0))^2)-1;
A2=2*pi*(f-1)*df*sqrt(m2*n2); %
r01=(sqrt(a0)-sqrt(ss.a1(k)*a0))/(sqrt(a0)+sqrt(ss.a1(k)*a0)); %反射系数
r12=(sqrt(ss.a1(k)*a0)-sqrt(ss.a2(k)*a0))/(sqrt(ss.a1(k)*a0)+sqrt(ss.a2(k)*a0));
r23=(sqrt(ss.a2(k)*a0)-sqrt(ss.a3(k)*a0))/(sqrt(ss.a2(k)*a0)+sqrt(ss.a3(k)*a0));
r1=r12*(1-r01^2)*exp(-2*A1*z1)*exp(-A0*z0); %广义反射系数
r2=r23*(1-r01^2)*(1-r12^2)*exp(-2*A1*z1)*exp(-2*A2*z2)*exp(-A0*z0);
re=(r1+r2)/2;ro=(r1-r2)/2; %偶奇分量
rf=ref_fft(f); % 反射系数频域序列
q1=2*re*cos(pi*(f-1)*df*ss.T(k))*cos(pi*(f-1)*df*(2*ss.t0(k)+ss.T(k)));
q2=2*ro*sin(pi*(f-1)*df*ss.T(k))*sin(pi*(f-1)*df*(2*ss.t0(k)+ss.T(k)));
q3=2*re*cos(pi*(f-1)*df*ss.T(k))*sin(pi*(f-1)*df*(2*ss.t0(k)+ss.T(k)));
q4=2*ro*sin(pi*(f-1)*df*ss.T(k))*cos(pi*(f-1)*df*(2*ss.t0(k)+ss.T(k)));
q=abs(real(rf)-q1-q2)+abs(imag(rf)+q3-q4);
Q2(k)=Q2(k)+q; %目标函数值(适应度函数)
% v1=c/sqrt(ss.a1(k));
% v2=c/sqrt(ss.a2(k));
% z0=0; %水平层厚度
% z1=v1*ss.t0(k)/2;
% z2=v2*ss.T(k)/2;
%
% A0=0; %衰减常数
% m1=(u0*ss.a1(k)*a0)/2;
% n1=sqrt(1+(ss.b1(k)/(2*pi*(f)*df*ss.a1(k)*a0))^2)-1;
% A1=2*pi*(f)*df*sqrt(m1*n1); %
% m2=(u0*ss.a2(k)*a0)/2;
% n2=sqrt(1+(ss.b2(k)/(2*pi*(f)*df*ss.a2(k)*a0))^2)-1;
% A2=2*pi*(f)*df*sqrt(m2*n2); %
%
% r01=(sqrt(a0)-sqrt(ss.a1(k)*a0))/(sqrt(a0)+sqrt(ss.a1(k)*a0)); %反射系数
% r12=(sqrt(ss.a1(k)*a0)-sqrt(ss.a2(k)*a0))/(sqrt(ss.a1(k)*a0)+sqrt(ss.a2(k)*a0));
% r23=(sqrt(ss.a2(k)*a0)-sqrt(ss.a3(k)*a0))/(sqrt(ss.a2(k)*a0)+sqrt(ss.a3(k)*a0));
%
% r1=r12*(1-r01^2)*exp(-2*A1*z1)*exp(-A0*z0); %广义反射系数
% r2=r23*(1-r01^2)*(1-r12^2)*exp(-2*A1*z1)*exp(-2*A2*z2)*exp(-A0*z0);
% re=(r1+r2)/2;ro=(r1-r2)/2; %偶奇分量
%
% rf=ref_fft(f); % 反射系数频域序列
% q1=2*re*cos(pi*(f)*df*ss.T(k))*cos(pi*(f)*df*(2*ss.t0(k)+ss.T(k)));
% q2=2*ro*sin(pi*(f)*df*ss.T(k))*sin(pi*(f)*df*(2*ss.t0(k)+ss.T(k)));
% q3=2*re*cos(pi*(f)*df*ss.T(k))*sin(pi*(f)*df*(2*ss.t0(k)+ss.T(k)));
% q4=2*ro*sin(pi*(f)*df*ss.T(k))*cos(pi*(f)*df*(2*ss.t0(k)+ss.T(k)));
% q=abs(real(rf)-q1-q2)+abs(imag(rf)+q3-q4);
% Q2(k)=Q2(k)+q; %目标函数值(适应度函数
end
end
%s.a1=9; %a,b根据实际情况要调整
%s.a2=25;
%s.a3=9;
%s.b1=0.02;
%s.b2=0.01;
%s.t0=7.8*10^-9;
%s.T=1*10^-9;
%Q1=10000000000000; %细细体会妙处(其实是偷懒不想编子函数)
n=0;
Q5=ones(1,364);
Q3=1e9;
Q4=1e9;
while T>1 %这是最外层的循环(温度),暂时不编写子函数
T0=T*r;
T=T0;
for i=1:M
%下为metropolis准则
for
评论16