function [Q1bar Q2bar]= CalculateQbar(Q, P31, P32, Pbest)

    PbestI=inv(Pbest);       
        Q1bar518=Q+Q*PbestI*P31+P31*PbestI*Q+Q*PbestI*Q+Q*PbestI*P31*PbestI*Q;
        Q2bar518=Q+Q*PbestI*P32+P32*PbestI*Q+Q*PbestI*Q+Q*PbestI*P32*PbestI*Q;
        Q1bar520=Q+Q*PbestI*P31+P31*PbestI*Q;
        Q2bar520=Q+Q*PbestI*P32+P32*PbestI*Q;
        Q1bar=Q+Q1bar518;
        Q2bar=Q+Q2bar518;
 
% Q1=4*Q;
% Q2=4*Q;

end