function MO_Animate(varargin)% This function generates objective space images showing why% sum-weighted optimizer can not find all non-dominated% solutions for non convex objective spaces in multi-ojective% optimization%% Guillaume JACQUENOTif nargin == 0 % Simu = 'Convex'; Simu = 'NonConvex'; save_pictures = true; interpreter = 'none';endswitch Simu case 'NonConvex' a = 0.1; b = 3; stepX = 1/200; stepY = 1/200; case 'Convex' a = 0.2; b = 1; stepX = 1/200; stepY = 1/200;end[X,Y] = meshgrid( 0:stepX:1,-2:stepY:2);F1 = X;F2 = 1+Y.^2-X-a*sin(b*pi*X);figure;grid on;hold on;box on;axis square;set(gca,'xtick',0:0.2:1);set(gca,'ytick',0:0.2:1);Ttr = get(gca,'XTickLabel');Ttr(1,:)='0.0';Ttr(end,:)='1.0';set(gca,'XTickLabel',[repmat(' ',size(Ttr,1),1) Ttr]);Ttr = get(gca,'YTickLabel');Ttr(1,:)='0.0';Ttr(end,:)='1.0';set(gca,'YTickLabel',[repmat(' ',size(Ttr,1),1) Ttr]);if strcmp(interpreter,'none') xlabel('f1','Interpreter','none'); ylabel('f2','Interpreter','none','rotation',0);else xlabel('f_1','Interpreter','Tex'); ylabel('f_2','Interpreter','Tex','rotation',0);endset(gcf,'Units','centimeters')set(gcf,'OuterPosition',[3 3 3+6 3+6])set(gcf,'PaperPositionMode','auto')[minF2,minF2_index] = min(F2);minF2_index = minF2_index + (0:numel(minF2_index)-1)*size(X,1);O1 = F1(minF2_index)';O2 = minF2';[pF,Pareto]=prtp([O1,O2]);fill([O1( Pareto);1],[O2( Pareto);1],repmat(0.95,1,3));text(0.45,0.75,'Objective space');text(0.1,0.9,'\leftarrow Optimal Pareto front','Interpreter','TeX');plot(O1( Pareto),O2( Pareto),'k-','LineWidth',2);plot(O1(~Pareto),O2(~Pareto),'.','color',[1 1 1]*0.8);V1 = O1( Pareto); V1 = V1(end:-1:1);V2 = O2( Pareto); V2 = V2(end:-1:1);O1P = O1( Pareto);O2P = O2( Pareto);O1PC = [O1P;max(O1P)];O2PC = [O2P;max(O2P)];ConvH = convhull(O1PC,O2PC);ConvH(ConvH==numel(O2PC))=[];c = setdiff(1:numel(O1P), ConvH);% Non convexO1PNC = O1PC(c);[temp, I1] = min(O1PNC);[temp, I2] = max(O1PNC);if ~isempty(I1) && ~isempty(I2) plot(O1PC(c),O2PC(c),'-','color',[1 1 1]*0.7,'LineWidth',2);endp1 = (V2(1)-V2(2))/(V1(1)-V1(2));hp = plot([0 1],[p1*(-V1(1))+V2(1) p1*(1-V1(1))+V2(1)]);delete(hp);Histo_X = [];Histo_Y = [];coeff = 0.02;Sq1 = coeff *[0 1 1 0 0;0 0 1 1 0];compt = 1;for i = 2:1:length(V1)-1 if ismember(i,ConvH) p1 = (V2(i+1)-V2(i-1))/(V1(i+1)-V1(i-1)); x_inter = 1/(1+p1^2)*(p1^2*V1(i)-p1*V2(i)); hp1 = plot([0 1],[p1*(-V1(i))+V2(i) p1*(1-V1(i))+V2(i)],'k'); % hp2 = plot([x_inter],[-x_inter/p1],'k','Marker','.','MarkerSize',8) hp3 = plot([0 x_inter],[0 -x_inter/p1],'k-'); hp4 = plot([x_inter 1],[-x_inter/p1 -1/p1],'k--'); hp5 = plot(V1(i),V2(i),'ko','MarkerSize',10); % Plot the square for perpendicular lines alpha = atan(-1/p1); Mrot = [cos(alpha) -sin(alpha);sin(alpha) cos(alpha)]; Sq_plot = repmat([x_inter;-x_inter/p1],1,5) + Mrot * Sq1; hp7 = plot(Sq_plot(1,:),Sq_plot(2,:),'k-'); Histo_X = [Histo_X V1(i)]; Histo_Y = [Histo_Y V2(i)]; hp6 = plot(Histo_X,Histo_Y,'k.','MarkerSize',10); w1 = p1/(p1-1); w2 = 1-w1; Fweight_sum = V1(i)*w1+w2*V2(i); Fweight_sum = floor(1e3*Fweight_sum )/1e3; w1 = floor(1000*w1)/1e3; str1 = sprintf('%.3f',w1); str2 = sprintf('%.3f',1-w1); str3 = sprintf('%.3f',Fweight_sum); if (strcmp(str1,'0.500')||strcmp(str1,'0,500')) && strcmp(Simu,'NonConvex') disp('Two solutions'); end title(['\omega_1 = ' str1 ' & \omega_2 = ' str2 ' & F = ' str3],'Interpreter','TeX'); axis([0 1 0 1]); file = ['Frame' num2str(1000+compt)]; if save_pictures saveas(gcf, file, 'epsc'); end compt = compt +1; pause(0.001); delete(hp1); delete(hp3); delete(hp4); delete(hp5); delete(hp6); delete(hp7); endenddisp(['Number of frames :' num2str(length(V1))]);return;function [A varargout]=prtp(B)% Let Fi(X), i=1...n, are objective functions% for minimization.% A point X* is said to be Pareto optimal one% if there is no X such that Fi(X)<=Fi(X*) for% all i=1...n, with at least one strict inequality.% A=prtp(B),% B - m x n input matrix: B=% [F1(X1) F2(X1) ... Fn(X1);% F1(X2) F2(X2) ... Fn(X2);% .......................% F1(Xm) F2(Xm) ... Fn(Xm)]% A - an output matrix with rows which are Pareto% points (rows) of input matrix B.% [A,b]=prtp(B). b is a vector which contains serial% numbers of matrix B Pareto points (rows).% Example.% B=[0 1 2; 1 2 3; 3 2 1; 4 0 2; 2 2 1;...% 1 1 2; 2 1 1; 0 2 2];% [A b]=prtp(B)% A =% 0 1 2% 4 0 2% 2 2 1% b =% 1 4 7A=[]; varargout{1}=[];sz1=size(B,1);jj=0; kk(sz1)=0;c(sz1,size(B,2))=0;bb=c;for k=1:sz1 j=0; ak=B(k,:); for i=1:sz1 if i~=k j=j+1; bb(j,:)=ak-B(i,:); end end if any(bb(1:j,:)'<0) jj=jj+1; c(jj,:)=ak; kk(jj)=k; endendif jj A=c(1:jj,:); varargout{1}=kk(1:jj);else warning([mfilename ':w0'],... 'There are no Pareto points. The result is an empty matrix.')endreturn;