How to create pareto front for bridge optimization program?

Hi guys,

I have problem with creating pareto front for this bridge optimization program. Below is the content of the task I am trying to do. As of now, I've been able to create a set of N random bridges with random positions of points P1 and P2 but honestly I don't know to apply Pareto front function to the rest of the code, so I could get this red line with most optimal bridge construction. Can you suggest some code I could use? Below you can see my code.

TASK

CODE:

s = 5; % width of the bridge
h = 10; % height of the bridge
F = 1000; % force 
pr = 230; % angle
N = 300; % number of variations

p = pr*2*pi/360; %conversion of degrees to radians

for i=0:N


% Coordinate points: 
P1x = rand*(s/2);
P1y = rand*h;
P2x = s/2 + rand*(s/2);
P2y = rand*h;

% Calculation of the length of the truss elements
d1 = s/2;
d2 = s/2;
d3 = sqrt(P1x.^2 + P1y.^2);
d4 = sqrt((s/2 - P1x).^2 + P1y.^2);
d5 = sqrt((P2x - s/2).^2 + P2y.^2);
d6 = sqrt((s - P2x).^2 + P2y.^2);
d7 = sqrt((P2x - P1x).^2 + (P1y - P2y).^2);

% Total length of elements (proportional to the cost of the structure)
d = s + d3 + d4 + d5 + d6 + d7;

% Calculation of angles between elements
k1 = acos( (d1^2 + d3^2 - d4^2) / (2*d1*d3) )*360/(2*pi);
k2 = acos( (d1^2 + d4^2 - d3^2) / (2*d1*d4) )*360/(2*pi);
k3 = acos( (d4^2 + d5^2 - d7^2) / (2*d4*d5) )*360/(2*pi);
k4 = acos( (d2^2 + d5^2 - d6^2) / (2*d2*d5) )*360/(2*pi);
k5 = acos( (d2^2 + d6^2 - d5^2) / (2*d2*d6) )*360/(2*pi);
k6 = acos( (d5^2 + d6^2 - d2^2) / (2*d5*d6) )*360/(2*pi);
k7 = acos( (d5^2 + d7^2 - d4^2) / (2*d5*d7) )*360/(2*pi);
k8 = acos( (d4^2 + d7^2 - d5^2) / (2*d4*d7) )*360/(2*pi);
k9 = acos( (d3^2 + d4^2 - d1^2) / (2*d3*d4) )*360/(2*pi);

% Sum of forces on the X axis of the entire system:
% EFx = cos(p)*F + RBx + RAx = 0 
% RBx + RAx = cos(p)*F
% RBx = RAx
RAx = F*cos(p) / 2;
RBx = RAx;


% Sum of moments relative to fixed bearing B:
% EMB = RAy*s - (s/2)*sin(p)*F = 0 
RAy = sin(p)*F/2;


% Sum of forces on the Y axis of the entire system:
% EFy = RAy + RBy - F*sin(p) = 0
RBy = F*sin(p) - RAy;



% Node 1.

% EFy = RAy + S15*sin(k1) = 0
S15 = -RAy / sin(k1);

% EFx = RAx + S12 + S15*cos(k1) = 0
S12 = -RAx - S15*cos(k1);



% Node 5.

% EFy = -sin(k1)*S15 - S52*sin(pi/2 - k1 - k9) - sin(pi - k9 -k8 -k1)*S54 =0
% EFx = -S15*cos(k1) + S54*cos(pi - k9 - k8 -k1) + S52*cos(pi - k1 - k9) =0
A = [sin(pi/2-k1-k9),sin(pi-k9-k8-k1);cos(pi-k1-k9),cos(pi-k9-k8-k1)];
B = [(-S15*sin(k1));S15*cos(k1)];
X = inv(A)*B;
S52 = X(1);
S54 = X(2);


% Node 2.

% EFy = S52*sin(k2) + S24*sin(k4) - F*sin(p)
S24 = ( F*sin(p) - S52*sin(k2) ) / ( sin(k4) );

% EFx = S23 - S12 - S52*cos(k2) + S24*cos(k4) + F*cos(p)
S23 = S12 + S52*cos(k2) - S24*cos(k4) - F*cos(p);


% Node 4.

% EFx = S43*sin(90-k5) - S54*cos(k5 - (pi -k6 -k7)) - S24*cos(pi-k6-k5) =0
S43 = ( S54*cos(k5 - (pi -k6 -k7)) + S24*cos(pi-k6-k5) ) / (sin(pi/2 -k5));

W = [abs(S12); abs(S15); abs(S52); abs(S54); abs(S24); abs(S23); abs(S43)];

Fmax = max(W)

 plot(d,Fmax,'-o')
 xlabel('dasda')
 ylabel('dsadada')
 axis([0 50 0 12000]); grid on
 hold on

end

[Widz Nieznany]

I've managed to find a code on matlab forum which I think is right for this problem but I don't know how to apply it. Maybe you guys have some ideas?

https://ch.mathworks.com/matlabcentral/fileexchange/50477-pareto-filtering?fbclid=IwAR3vTtdQTs6Ogy0Vlp5DtQ1nU3x_GM53x1mhmN03kRiTUmGa5X-hyqSwd1E

function [ p, idxs] = paretoFront( p )
% Filters a set of points P according to Pareto dominance, i.e., points
% that are dominated (both weakly and strongly) are filtered.
%
% Inputs: 
% - P    : N-by-D matrix, where N is the number of points and D is the 
%          number of elements (objectives) of each point.
%
% Outputs:
% - P    : Pareto-filtered P
% - idxs : indices of the non-dominated solutions
%
% Example:
% p = [1 1 1; 2 0 1; 2 -1 1; 1, 1, 0];
% [f, idxs] = paretoFront(p)
%     f = [1 1 1; 2 0 1]
%     idxs = [1; 2]

[i, dim] = size(p);
idxs = [1 : i]';
while i >= 1
    old_size = size(p,1);
    indices = sum( bsxfun( @ge, p(i,:), p ), 2 ) == dim;
    indices(i) = false;
    p(indices,:) = [];
    idxs(indices) = [];
    i = i - 1 - (old_size - size(p,1)) + sum(indices(i:end));
end    

end

Answer 1

[Francisco Martos Barrachina]

Hi there!

I imagine it might be to late for you, but it is never too late for an answer!

What you first need is a matrix (let's call it ov for objective value) comprising the values of the two coordinates you have in that 2D plot with highest tension and price. The two coordinates are your objectives and every point is a viable (feasible) solution. It seems you want to minimize both, which is fine.

That ov matrix should look like this:

% Random values for ov
% The columns represent the two objectives
% Column one could be Price (x axys in your plot)
% Column two colud be Tension (y axys in your plot)
% The rows represent every bridge evaluated
ov= [1,3; % bridge 1 has a price of 1 and a tension of 3
     2,5; % ...
     4,1; % ... 
     3,2; % ...
     5,1; % ...
     3,2] % brige 6 has a price of 3 and a tension of 2

With that, you can use it as input in a pareto function. I have a custom one that I will happily share. In mine, all objectives are supposed to be minimized (usual in multiobjective optimization), so if you have an objective to maximize, you can just change it from positive to negative. My function is actually two functions nested inside each other. Variable names are in Spanish. You can find them here https://github.com/curromb20/matlab_functions.git.

matrizObjetivos, which in our little example would be 'ov' is the only input, again with points in rows and objectives in columns

function [matrizEficientes,indiceEficientes] = paretoGlobal_min(matrizObjetivos)
% Inputs are:
% matrizObjetivos -> a matrix of objective values 
% The number of rows is the number of feasible solutions
% The number of columns is the number of objectives
% Outputs are:
% matrizEficientes -> objective values of the efficient solutions
% indiceEficientes -> index of the efficient solutions

matrizEficientes=[];
indiceEficientes=[];
[numSol,~]=size(matrizObjetivos);

for i=1:numSol
    
    % here we use esPareto_min 
    [vectorEficientes,logico] = esPareto_min(matrizEficientes,matrizObjetivos(i,:));
    
    
    if isempty(vectorEficientes)
        matrizEficientes=matrizObjetivos(i,:);
        indiceEficientes=i;
    elseif logico
        matrizEficientes=matrizEficientes(logical(vectorEficientes),:);
        indiceEficientes=indiceEficientes(logical(vectorEficientes));
        matrizEficientes(end+1,:)=matrizObjetivos(i,:);
        indiceEficientes(end+1)=i;
    end
    
end
end

The function we additionally use is esPareto_min, that evaluates if any new solution (compared with a set of them) is efficient.

function [vectorIndEficientes,logico] = esPareto_min(matrizEficientes,nuevaSol)
% Inputs are:
% matrizEficientes -> which is a set of the objective values of known efficient solutions
% nuevaSol -> the objective values of a new solution
% They both have as many columns as there are objectives
% We want to check if the new solution is efficient
% If it is, we check if it dominates any of the current set of efficients
% If it does, the newly dominated solutions get thrown out of matrizEficientes 
% Every objective is to be minimized
% Outputs:
% vectorIndEficientes -> vector with index of solutions that are still Efficient
% logico -> true if nuevaSol is efficient

logico=true;

if isempty(matrizEficientes)
    vectorIndEficientes = [];
else
    % First check if it coincides exactly with an efficient solution
    
    [numSolEfic,numObjetivos]=size(matrizEficientes);
    vectorIndEficientes = ones(numSolEfic,1); %1 = efficient.
    
    nuevaIgual=matrizEficientes==nuevaSol;
    % to each solution in each objective
    
    nuevaIgualGlobal=(sum(nuevaIgual,2)==numObjetivos);
    % to each solution in every objective
    
    if nuevaIgualGlobal>0
        % if it is equal to an efficient solution
        % just leave the vector of efficient solutions as it is
        % and have logico as true
    else
        % see if it dominates 
        nuevaDomina=matrizEficientes>=nuevaSol;
        % to each solution in each objective
        
        nuevaDominaGlobal=(sum(nuevaDomina,2)==numObjetivos);
        % to each solution in every objective
        
        % see if it is dominated
        nuevaEsDominada=nuevaSol>=matrizEficientes;
        % by each solution in each objective
        
        nuevaEsDominadaGlobal=(sum(nuevaEsDominada,2)==numObjetivos);
        % by each solution in every objective
        
        if sum(nuevaEsDominadaGlobal)>0
            % if it is dominated
            logico=false;
        elseif sum(nuevaDominaGlobal)>0
            % if it dominates
            vectorIndEficientes=1-nuevaDominaGlobal;
            % gives back indexes of solutions that are still efficient
        end
    end
end
end

Using this should be pretty straightforward. And I don't think the plot is a problem. One tip though, if plotting is problematic, is sorting the pareto efficient solutions by one of the objectives (as you have two) and it will be reversely sorted by the other, so you can link the resulting dots by this order.

I hope this can be of help!