How to maximize determinant of a matrix?

15 visualizaciones (últimos 30 días)
Ammy
Ammy el 4 de Ag. de 2021
Comentada: Walter Roberson el 4 de Ag. de 2021
I want to apply ABC OPTIMIZATION algorithm to maximize determinant of a matrix.
I have my matrix in the form
X=A(nxn) +k.*B(nxn)
where A and B are nxn matrices, and k is any number between 0 to 10.
I want to optimize k with ABC algorithm such that det(X) is maximum, in other words I have to find such k for which det (X) is maximum.
The code for ABC is attached I need help regarding to define objective function for this.
%
clc;
clear;
close all;
%% Problem Definition
CostFunction=@(x) ObjectiveFunction(x); % Cost Function
nVar=1; % Number of Decision Variables
VarSize=[1 nVar]; % Decision Variables Matrix Size
VarMin=0; % Decision Variables Lower Bound
VarMax= 10; % Decision Variables Upper Bound
%% ABC Settings
MaxIt=200; % Maximum Number of Iterations
nPop=100; % Population Size (Colony Size)
nOnlooker=nPop; % Number of Onlooker Bees
L=round(0.6*nVar*nPop); % Abandonment Limit Parameter (Trial Limit)
a=1; % Acceleration Coefficient Upper Bound
%% Initialization
% Empty Bee Structure
empty_bee.Position=[];
empty_bee.Cost=[];
% Initialize Population Array
pop=repmat(empty_bee,nPop,1);
% Initialize Best Solution Ever Found
BestSol.Cost=inf;
% Create Initial Population
for i=1:nPop
pop(i).Position=unifrnd(VarMin,VarMax,VarSize);
pop(i).Cost=CostFunction(pop(i).Position);
if pop(i).Cost<=BestSol.Cost
BestSol=pop(i);
end
end
% Abandonment Counter
C=zeros(nPop,1);
% Array to Hold Best Cost Values
BestCost=zeros(MaxIt,1);
%% ABC Main Loop
for it=1:MaxIt
% Recruited Bees
for i=1:nPop
% Choose k randomly, not equal to i
K=[1:i-1 i+1:nPop];
k=K(randi([1 numel(K)]));
% Define Acceleration Coeff.
phi=a*unifrnd(-1,+1,VarSize);
% New Bee Position
newbee.Position=pop(i).Position+phi.*(pop(i).Position-pop(k).Position);
% Evaluation
newbee.Cost=CostFunction(newbee.Position);
% Comparision
if newbee.Cost<=pop(i).Cost
pop(i)=newbee;
else
C(i)=C(i)+1;
end
end
% Calculate Fitness Values and Selection Probabilities
F=zeros(nPop,1);
MeanCost = mean([pop.Cost]);
for i=1:nPop
F(i) = exp(-pop(i).Cost/MeanCost); % Convert Cost to Fitness
end
P=F/sum(F);
% Onlooker Bees
for m=1:nOnlooker
% Select Source Site
i=RouletteWheelSelection(P);
% Choose k randomly, not equal to i
K=[1:i-1 i+1:nPop];
k=K(randi([1 numel(K)]));
% Define Acceleration Coeff.
phi=a*unifrnd(-1,+1,VarSize);
% New Bee Position
newbee.Position=pop(i).Position+phi.*(pop(i).Position-pop(k).Position);
% Evaluation
newbee.Cost=CostFunction(newbee.Position);
% Comparision
if newbee.Cost<=pop(i).Cost
pop(i)=newbee;
else
C(i)=C(i)+1;
end
end
% Scout Bees
for i=1:nPop
if C(i)>=L
pop(i).Position=unifrnd(VarMin,VarMax,VarSize);
pop(i).Cost=CostFunction(pop(i).Position);
C(i)=0;
end
end
% Update Best Solution Ever Found
for i=1:nPop
if pop(i).Cost<=BestSol.Cost
BestSol=pop(i);
end
end
% Store Best Cost Ever Found
BestCost(it)=BestSol.Cost;
% Display Iteration Information
disp(['Iteration ' num2str(it) ': Best Cost = ' num2str(BestCost(it))]);
end
%% Results
figure;
%plot(BestCost,'LineWidth',2);
semilogy(BestCost,'LineWidth',2);
xlabel('Iteration');
ylabel('Best Cost');
grid on;

Respuesta aceptada

Walter Roberson
Walter Roberson el 4 de Ag. de 2021
CostFunction = @(x) -det(x); % Cost Function
  4 comentarios
Ammy
Ammy el 4 de Ag. de 2021
Thank you very much.
Walter Roberson
Walter Roberson el 4 de Ag. de 2021
Note that for the system A+k*B that the optimal K can be calculated without searching, if you have the symbolic toolbox
syms k
C = A * k.*B;
dC = diff(det(C),k);
k_candidates = solve(dC, k);
LB = 0; UB = 10;
mask = imag(k_candidates) == 0 & real(k_candidates) >= LB & real(k_candidates) <= UB;
selected_k = [k_candidates(mask), LB, UB]
selected_det = arrayfun(@(K) det(subs(C, k, K))), selected_k)
[best_det, idx] = max(selected_det);
best_k = selected_k(idx)
vpa(best_det)
vpa(best_k)

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Get Started with Optimization Toolbox en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by