Path: news.mathworks.com!not-for-mail
From: "kuffal al" <kuff_logg@yahoo.com.my>
Newsgroups: comp.soft-sys.matlab
Subject: Re: optimization of fiber bragg grating using nelder-mead simplex
Date: Sun, 4 Apr 2010 16:43:06 +0000 (UTC)
Organization: mmu
Lines: 273
Message-ID: <hpafiq$s2m$1@fred.mathworks.com>
References: <hpa96d$6t1$1@fred.mathworks.com> <hpaec3$e6o$1@fred.mathworks.com>
Reply-To: "kuffal al" <kuff_logg@yahoo.com.my>
NNTP-Posting-Host: webapp-02-blr.mathworks.com
Content-Type: text/plain; charset=UTF-8; format=flowed
Content-Transfer-Encoding: 8bit
X-Trace: fred.mathworks.com 1270399386 28758 172.30.248.37 (4 Apr 2010 16:43:06 GMT)
X-Complaints-To: news@mathworks.com
NNTP-Posting-Date: Sun, 4 Apr 2010 16:43:06 +0000 (UTC)
X-Newsreader: MATLAB Central Newsreader 2298588
Xref: news.mathworks.com comp.soft-sys.matlab:623364

i actually dont really understand about the objective function bcoz there is a lots of equation. so which one should i choose or how should i create an objective function.


> 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=@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))) 
> >