Prev: counting backwards
Next: Remapping Matrix Values (3d to 2d) for Convolution
From: Alan Weiss on 4 Apr 2010 12:22
You might want to use some built-in MATLAB functions.
fminsearch performs the Nelder-Mead simplex algorithm.
It accepts a single vector of inputs, say x = [x(1),x(2),x(3)], and
works on your objective function that should take x and return f(x),
where f(x) is the function to minimize. For your problem, take x =
[L,A,dn], and write your objective function in terms of x.

It seems to me that your code is doing some very inefficient
root-finding as well. Take a look at the fzero function for a solid
root-finding implementation.

fminsearch reference:
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fminsearch.html

fminsearch example:
http://www.mathworks.com/access/helpdesk/help/techdoc/math/bsgpq6p-10.html

fzero reference:
http://www.mathworks.com/access/helpdesk/help/techdoc/ref/fzero.html

fzero example:
http://www.mathworks.com/access/helpdesk/help/techdoc/math/bsgpq6q-44.html

Alan Weiss
MATLAB mathematical toolbox documentation

kuffal al wrote:
> 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: counting backwards
Next: Remapping Matrix Values (3d to 2d) for Convolution