How to maximize determinant of a matrix?

15 vues (au cours des 30 derniers jours)
Ammy
Ammy le 4 Août 2021
Commenté : Walter Roberson le 4 Août 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;

Réponse acceptée

Walter Roberson
Walter Roberson le 4 Août 2021
CostFunction = @(x) -det(x); % Cost Function
  4 commentaires
Ammy
Ammy le 4 Août 2021
Thank you very much.
Walter Roberson
Walter Roberson le 4 Août 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)

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Get Started with Optimization Toolbox dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by