Prev: Remapping Matrix Values (3d to 2d) for Convolution
Next: Watershed algorithm
From: kuffal al on 4 Apr 2010 12:43
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=(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)))
> >
From: James Allison on 6 Apr 2010 16:16
It seems that this might be formulated more naturally as a constrained
design optimization problem, i.e., you need to come up with both an
objective function AND constraint functions that fit the following form:

minimize f(x)
subject to g(x) <= 0
with respect to x

where x is the optimization vector (the independent variables you are
trying to adjust), and g may be a vector-valued function. What should be
the objective and what should be the constraints requires some
domain-specific knowledge about your application. What is most important
about the design? Are you trying to minimize or maximize something
(i.e., and objective)? Are there some limits that should not be exceeded
(i.e., constraints)?

fminsearch (Nelder-Mead) only applies to unconstrained problems, so it
won't work if you formulate your problem as a constrained optimization
problem. fmincon works well if the functions are smooth. patternsearch
would be the next thing to try if the functions are not smooth. The
genetic algorithm (ga) can solve very difficult optimization problems,
but requires a large number of function evaluations.

-James

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