MATLAB Examples

## Contents

```% Set up problem % % % domain x_min=[-1;-1.00]; %min values for each dimension x_max=[ 1; 1.25]; %max values for each dimension % % % functions, label one Y and Y_Grad, and another set J and J_Grad % straight up parabola J =@(x) +((x(1)-0).^2+(x(2)-0.2).^2); J_Grad=@(x) [2*(x(1)-0);2*(x(2)-0.2)]; % linear in x1 and quadratic in x2 % Y=@(x) 1.25*(x(1)-0.2).*(x(2).^2); % Y_Grad=@(x) [1.25*(x(2).^2);1.25*(x(1)-0.2).*2*x(2)]; % linear in 1d and parabola in other dimension % Y =@(x) (x(1).^1).*(sin(2*pi/2*x(2))); % Y_Grad=@(x) [(sin(2*pi/2*x(2)));2*pi/2*(x(1).^1).*(cos(2*pi/2*x(2)))]; % Multimodal gausian, peaks(a1*2,a2*2)/7, equation is used to show gradient Y =@(x) (3*(1-2*x(1))^2*exp(-4*x(1)^2-(2*x(2)+1)^2)+... -exp(-(2*x(1)+1)^2-4*x(2)^2)/3+... -10*(2*x(1)/5-8*x(1)^3-32*x(2)^5)*exp(-4*x(1)^2-4*x(2)^2))/7; Y_Grad=@(x) [-(12+24*x(1)*(1-2*x(1)))*(1-2*x(1))*exp(-4*x(1)^2-(2*x(2)+1)^2)+... 4*(2*x(1)+1)*exp(-(2*x(1)+1)^2-4*x(2)^2)/3+... 10*(-2/5+24*x(1)^2+8*x(1)*(2*x(1)/5-8*x(1)^3-32*x(2)^5))*exp(-4*x(1)^2-4*x(2)^2);... 3*(1-2*x(1))^2*(-4*(2*x(2)+1))*exp(-4*x(1)^2-(2*x(2)+1)^2)+... +8*x(2)*exp(-(2*x(1)+1)^2-4*x(2)^2)/3+... +10*(160*x(2)^4+8*x(2)*(2*x(1)/5-8*x(1)^3-32*x(2)^5))*exp(-4*x(1)^2-4*x(2)^2)]/7; ```

## optimization problem statement and parameters

```% ProbStrut.J=@(x) x(1).^2+x(2).^2; % cost function % ProbStrut.J_Grad=@(x) [2*x(1);2*x(2)]; % gradient of the cost function % ProbStrut.y=@(x) x(2); % output function % ProbStrut.y_Grad=@(x) [0;1]; % gradient of the cost function % ProbStrut.y=@(x) x(1)+x(2); % output function % ProbStrut.y_Grad=@(x) [1;1]; % gradient of the cost function ProbStrut.J =J; ProbStrut.y =Y; ProbStrut.J_Grad=J_Grad; ProbStrut.y_Grad=Y_Grad; ProbStrut.sat = @(x,beta) x/max(norm(x),beta); %saturation function used % domain, used to size initial population ProbStrut.x_min=x_min; ProbStrut.x_max=x_max; % step parameters ProbStrut.beta_J=0.05; % saturation length for cost gradient ProbStrut.beta_y=0.025; % saturation length for output gradient ProbStrut.beta_Delta=0.5; % saturation length for step direction ProbStrut.beta_delta=0.88; % step length reduction factor for Armijo condition ProbStrut.sigma=0.2; % Armijo rule factor ProbStrut.gamma_J=0.05; % max step length for cost reduction steps ProbStrut.gamma_y=0.03; % max step length for output change steps % cluster parameters ProbStrut.knots_min=7.5e-6; %If new knot is within this, cluster is terminated ProbStrut.knots_max=0.1;%If new knot is further than this, cluster is terminated ProbStrut.knots_itter_max=150; %If new knot does not converge within this, cluster is terminated % other parameters ProbStrut.pop_init_size=60; % number of initial agents to search space ProbStrut.Delta_min=1e-4; % if step length is less, agent is converged ProbStrut.pop_dist_kill=0.15; % points within this distance are terminated ProbStrut.pop_itter_max=400; % number of itterations for population steps ProbStrut.BoundViolateOK=1e-10; %so there isn't numerical issues with bounds ```

## set up function to watch particles move

plot function is optional argument first argument is a list of the population (columns are agents) second argument is a style for different phases: 'b.' for moving, 'r+' for cluster,'c*' for cluster start point and 'go' for agent removed due to proximity.

```% give the figure a background figure(61);clf; [temp_xx,temp_yy]=meshgrid(linspace(ProbStrut.x_min(1),ProbStrut.x_max(1),40),... linspace(ProbStrut.x_min(2),ProbStrut.x_max(2),55)); YY=zeros(size(temp_xx));JJ=zeros(size(temp_xx)); for ella=1:numel(temp_xx) YY(ella)=ProbStrut.y([temp_xx(ella);temp_yy(ella)]); JJ(ella)=ProbStrut.J([temp_xx(ella);temp_yy(ella)]); end contour(temp_xx,temp_yy,YY);hold on; %need hold to be on to see progression contour(temp_xx,temp_yy,JJ); PlotStep_axe=gca; % this ensures same figure is used each time title({'Progress of optimization';'Legend: b.: path, r+: cluster, ';... 'c*: Initial cluster point and go for killed path'}); ProbStrut.Plot_Step=@(Pops,Symb) plot(PlotStep_axe,Pops(1,:),Pops(2,:),Symb); ProbStrut.Plot_final=@(Pops,Symb) pause(0.1); %just to slow it down to watch ```

## Initialize population

```Pop_I=bsxfun(@plus,bsxfun(@times,rand(length(ProbStrut.x_min),ProbStrut.pop_init_size),... ProbStrut.x_max-ProbStrut.x_min),ProbStrut.x_min); ```

## Call optimization

```InvFuntions=InvFun(Pop_I,ProbStrut); % InvFuntions=InvFun([0;0],ProbStrut) ```

## Comparison of all inverse functions

```temp_color=jet(length(InvFuntions)); figure(41);clf; for carly=1:length(InvFuntions) for dani=1:length(ProbStrut.x_min) subplot(length(ProbStrut.x_min)+1,1,dani) hold on;grid on; plot(InvFuntions{carly}.y,InvFuntions{carly}.x(dani,:),'-+','color',temp_color(carly,:)); ylabel(['Input x_' num2str(dani)]); end subplot(length(ProbStrut.x_min)+1,1,length(ProbStrut.x_min)+1) hold on;grid on; plot(InvFuntions{carly}.y,InvFuntions{carly}.J,'-+','color',temp_color(carly,:)); ylabel('J, cost'); xlabel('y, output'); end ```

## Accuracy and optimality of all inverse functions

```Pnts_on=50; Pnts_off=200; Pnts_off_baseNoise=0.01; %additional noise for near points Pnts_mrk=5; %marker size for dots in plots for carly=1:length(InvFuntions) figure(100+carly);clf; % points on the inverse function Pnts_on_test_y=rand(1,Pnts_on)*(InvFuntions{carly}.y(end)-InvFuntions{carly}.y(1))+InvFuntions{carly}.y(1); Pnts_on_test_x=pchip(InvFuntions{carly}.y,InvFuntions{carly}.x,Pnts_on_test_y); Pnts_on_test_y_act=zeros(size(Pnts_on_test_y)); Pnts_on_test_J=zeros(size(Pnts_on_test_y)); for amy=1:size(Pnts_on_test_x,2) Pnts_on_test_y_act(amy)=ProbStrut.y(Pnts_on_test_x(:,amy)); Pnts_on_test_J(amy)=ProbStrut.J(Pnts_on_test_x(:,amy)); end % nearby points % using testing points should be more uniform in output... maybe temp_cov=cov(Pnts_on_test_x.').'+Pnts_off_baseNoise*eye(length(ProbStrut.x_min)); Pnts_off_test_x=bsxfun(@plus,sqrtm(temp_cov)*randn(length(ProbStrut.x_min),Pnts_off),... mean(Pnts_on_test_x,2)); Pnts_off_test_x=bsxfun(@max,bsxfun(@min,Pnts_off_test_x,ProbStrut.x_max),ProbStrut.x_min); Pnts_off_test_y=zeros(1,size(Pnts_off_test_x,2)); Pnts_off_test_J=zeros(1,size(Pnts_off_test_x,2)); for amy=1:size(Pnts_off_test_x,2) Pnts_off_test_y(amy)=ProbStrut.y(Pnts_off_test_x(:,amy)); Pnts_off_test_J(amy)=ProbStrut.J(Pnts_off_test_x(:,amy)); end subplot(1,2,1) plot(Pnts_off_test_x(1,:),Pnts_off_test_x(2,:),'g.','markersize',Pnts_mrk); hold on; plot(InvFuntions{carly}.x(1,:),InvFuntions{carly}.x(2,:),'r+'); plot(Pnts_on_test_x(1,:),Pnts_on_test_x(2,:),'m.','markersize',Pnts_mrk); contour(temp_xx,temp_yy,YY); contour(temp_xx,temp_yy,JJ); xlabel('x_1');ylabel('x_2'); title('Location of test and near points'); legend({'Near points','Inverse function','Test points'},'location','best') axis([ProbStrut.x_min(1),ProbStrut.x_max(1),ProbStrut.x_min(2),ProbStrut.x_max(2)]); subplot(2,2,2) plot([InvFuntions{carly}.y(end),InvFuntions{carly}.y(1)],[0,0],'r-');hold on; plot(Pnts_on_test_y,Pnts_on_test_y-Pnts_on_test_y_act,'m.','markersize',Pnts_mrk); xlabel('Output, y');ylabel('Error, predicted-actual'); title('Accuracy of test points'); subplot(2,2,4) plot(Pnts_off_test_y,Pnts_off_test_J,'g.','markersize',Pnts_mrk);hold on; plot(InvFuntions{carly}.y,InvFuntions{carly}.J,'r+'); plot(Pnts_on_test_y,Pnts_on_test_J,'m.','markersize',Pnts_mrk); xlabel('Output, y');ylabel('Cost'); title('Optimality of test points'); legend({'Near points','Inverse function','Test points'}) end ```