Prev: imagesc
Next: Resetting for loop
From: kuffal al on 4 Apr 2010 10:54
can someone help me to optimize the fiber bragg grating(fbg) using nelder mead simplex algoritm. i have both codes (both codes working well) the fbg n nelder mead but i fail to integrate them. the parameter that i want to optimize are the length (L), period (A) and the induced index change (dn). so here are the codes

(1)the fbg code
%==========================================================================
% matlab script for solving the coupled-mode equation using
% the transfer matrix method
%==========================================================================

clear all
close all
clc

%==========================================================================
% fibre simulation parameters
%==========================================================================

L = 0.0100; % Grating Length

dn = 1.0e-4; % Induced index change

design_lambda = 1550e-9; % design wavelength

%==========================================================================
% finding neff(effective reflective index) for the designed fibre
%==========================================================================

format long
n=40;
a=0.5;
increment=2;
tolerance=1e-4;

nco=1.447; % nco-core index
ncl=1.444; % ncl-cladding index
xstart=((ncl-0.0112)+nco)/2;
x=xstart;
dx=increment;

%==========================================================================
% to find the root of neff(effective reflective index)
%==========================================================================


for m=1:n
neff= sqrt((ncl^2)+(a*x)*(nco^2-ncl^2));
while dx/x>tolerance
if neff~=sqrt((ncl^2)+(a*(x+dx))*(nco^2-ncl^2));
dx=dx/2;
else
x=x+dx;
end
end
route(m)=x;
dx=increment;
x=1.0001*x;
end

neff=route(40:end);

V =(2*pi/design_lambda)*a*sqrt(nco^2-ncl^2); % normalize frequency

A=design_lambda/(2*neff); % Grating period 5.36e-7;

lambda=(1549.50e-9:0.005e-9:1550.50e-9); % Spectra range


%==========================================================================
% Coupled-mode Theory
%==========================================================================

v=1; % visibility, v assumed to be 1

k=(pi./lambda)*v*dn; % AC coupling coefficient

% DC(period averaged) coupling coefficient

% for stronger grating, sigma_dc assumed to be 0

sigma_dc =0; %(2*pi./lambda).*dn;

% general DC coupling coefficient

sigma = 2*pi*neff.*(1./lambda - 1/(2*neff*A)) + sigma_dc;


%==========================================================================
% Power Reflection coefficients(R)
%==========================================================================

K2 = k.*k;

sigma2 = sigma.*sigma;

difference = K2-sigma2;

up = sinh(L*sqrt(difference)).^2;

down = cosh(L*sqrt(difference)).^2 - sigma2./K2;

R = up./down;

Rt = transpose(R);

%==========================================================================
% saving the points
%==========================================================================

save 'fileToRead1.txt' Rt -ascii

%==========================================================================
% plots for reflection spectra for Bragg reflector
%==========================================================================

figure(1)

plot(lambda,R,'b');

axis( [1549.50e-9, 1550.50e-9, 0, 1]);

title('Reflection Spectra of Bragg grating');

xlabel('Wavelength(m)');

ylabel('Reflectivity(p.u)');

figure(2)
stairs(lambda,R,'r');figure(gcf)
axis( [1549.50e-9, 1550.50e-9, 0, 1]);

--------------------------------------------------------------------------------------------------------------------

(2) the nelder-mead code

%*******************************************
%* Nelder and Mead Technique *
%* Siamak Faridani *
%* University of Oklahoma *
%* 16 Oct 2005 *
%* www.PitchUp.com *
%*******************************************


clc;
clear all;
close all;
N=2;
alpha=2;
beta=.25;
gamma=2.5;
Precision=.00001

%Insert your starting point and the function here
x01=1.5505*1e-6;
x02=1.550e-6;
func=(a)Modified_fbg


%do not modify the following lines
%==============================

x=[[x01,x02]];



delta1=(((N+1)^.5+N-1)/(N*sqrt(2)))*alpha;

delta2=(((N+1)^.5-1)/(N*sqrt(2)))*alpha
x=[x;[x(1,1)+delta1,x(1,2)+delta2];[x(1,1)+delta2,x(1,2)+delta1]];
x_archive=[x [func(x(1,1),x(1,2));func(x(2,1),x(2,2));func(x(3,1),x(3,2))]];
l=0
x
disp('start================');
Precision_inverse=Precision^-1;
while ~(isequal(round(x(1,:)*Precision_inverse),round(x(2,:)*Precision_inverse),round(x(3,:)*Precision_inverse)))
l=l+1;
values=[func(x(1,1),x(1,2));func(x(2,1),x(2,2));func(x(3,1),x(3,2))];

sorted=[find(values==max(values)),6-find(values==max(values))-find(values==min(values)),find(values==min(values))];
xc=(1/N)*(x(sorted(1,2),:)+x(sorted(1,3),:));
disp('trying Normal Reflection');
xnew=x(sorted(1,1),:)+(1+alpha)*(xc-x(sorted(1,1),:))
isitok=0;
func(xnew(1,1),xnew(1,2))
while isitok==0

if ((func(xnew(1,1),xnew(1,2))>values(sorted(1,3),:)) && (func(xnew(1,1),xnew(1,2))<values(sorted(1,2),:)))
disp('f(l)<f(x_new)<f(g)')
isitok=1;
elseif (func(xnew(1,1),xnew(1,2))<values(sorted(1,3),:))
disp('lets Expand...');
xnew1=x(sorted(1,1),:)+(1+gamma)*(xc-x(sorted(1,1),:))
func(xnew1(1,1),xnew1(1,2))
if (func(xnew1(1,1),xnew1(1,2))<func(xnew(1,1),xnew(1,2)))
xnew=xnew1
end
isitok=1;
elseif (func(xnew(1,1),xnew(1,2))>values(sorted(1,2),:)) && (func(xnew(1,1),xnew(1,2))>values(sorted(1,1),:))
disp('lets Contract with teta=-beta...');
xnew1=x(sorted(1,1),:)+(1-beta)*(xc-x(sorted(1,1),:))
func(xnew1(1,1),xnew1(1,2))
if (func(xnew1(1,1),xnew1(1,2))<func(xnew(1,1),xnew(1,2)))
xnew=xnew1
end
isitok=1;
elseif (func(xnew(1,1),xnew(1,2))>values(sorted(1,2),:)) && (func(xnew(1,1),xnew(1,2))<values(sorted(1,1),:))
disp('lets Contract with teta=beta...');
xnew1=x(sorted(1,1),:)+(1+beta)*(xc-x(sorted(1,1),:))
func(xnew1(1,1),xnew1(1,2))
if (func(xnew1(1,1),xnew1(1,2))<func(xnew(1,1),xnew(1,2)))
xnew=xnew1
end
isitok=1;
end
end
x(sorted(1,1),:)=xnew;
x
disp('Finished================');
x_archive=[x_archive;[x [func(x(1,1),x(1,2));func(x(2,1),x(2,2));func(x(3,1),x(3,2))]]];
end

subplot(1,2,1)
plot(x_archive(:,1),x_archive(:,2),'o-r')
subplot(1,2,2)
patch(x_archive(:,1),x_archive(:,2),x_archive(:,3),'r')
view(3)

values=[func(x(1,1),x(1,2));func(x(2,1),x(2,2));func(x(3,1),x(3,2))];

sorted=[find(values==max(values)),6-find(values==max(values))-find(values==min(values)),find(values==min(values))];
disp(sprintf('Number of iteration= %10d\n',l))
disp(sprintf('minimum of f(x) is %10.5f that occurs at x=[%10.5f,%10.5f]\n',min(values),x(find(values==min(values)),1),x(find(values==min(values)),2)))
 | 
Pages: 1
Prev: imagesc
Next: Resetting for loop