Using the given code, how can I optimize a matrix to minimize a cost function?

Given S1, (1,K) vector, I want to optimize B (N,M) matrix to minimize the following cost function:

Subject to:

Where:

S2 is (1,K) vector and a function of matrix B.

S2 can be calculated after optimizing matrix B using the following system parameters and equations:

clc;
clear;

% Given system parameters:
N = 2;
K = 4;
M=2;
C_l = 4;
H = [0.1185    0.2811; 0.3550    0.8224; 0.3260    0.9644; 0.5333    0.6083]; % 4*2 matrix
A = [-2 1; -1 1]; % 2*2 matrix
C = [7 -3; 7 -3; -2 1; -2 1]; % 4*2 matrix
P = [25000000   0; 0  25000000]; % 4*4 matrix
S1 = [3.1683    3.1686    1.8716    1.8898]; % 1*4 vector
S2 = zeros(1,K); % intial value
B = zeros(N,M); % intial value

% How can we  optimize the value of B matrix to achieve our goal? 

%calculate S2 from B and the other given inputs
for j=1:1:N
    d(j) = (B(j,:)*P*B(j,:)')/((2^(2*C_l))-(norm(A(:,j))^2));
end
D_d = diag(d);

for i=1:1:K   
    V_d(i)=C(i,:)*P*B'*H(i,:)'*inv(1+H(i,:)*(A'*D_d*A+B*P*B')*H(i,:)');
    sigma_d(i)=norm((V_d(i)*H(i,:)*B-C(i,:))*(P^(1/2)))^2+(V_d(i)^2)*(1+H(i,:)*A'*D_d*A*H(i,:)');
    
    S2(i)=0.5*log2((P(1,1))/sigma_d(:,i));   
end

[Stephan]

Do you have at least one feasible B, that satisfies the constraints?

[Adi Nor]

No

[Rik]

Can you at least write the cost function?

[Adi Nor]

The cost function is the first equation in my question

because I want each element in vector S2 to be equal that in the given S1. for example: I want to make S2(1,1) = S1 (1,1) and so on till S2(1,K) = S1(1,K). So, I think the above cost function may be written in MATLAB as:

sum(S2 - S1).^2

But, I am not sure.

[Rik]

I edited your code. I implemented your constraint inside the cost function.

I am using fminsearch, which can be sensitive to a local minimum depending on the initial guess. In this case that doesn't matter, because it can't find any valid solution. Are you sure it exists? Are you sure you did not make any mistake when implementing the mathematics in Matlab? Either way, the code below could be helpful.

clc;clear;

% Given system parameters:
N = 2;
K = 4;
M=2;
C_l = 4;
H = [0.1185    0.2811; 0.3550    0.8224; 0.3260    0.9644; 0.5333    0.6083]; % 4*2 matrix
A = [-2 1; -1 1]; % 2*2 matrix
C = [7 -3; 7 -3; -2 1; -2 1]; % 4*2 matrix
P = [25000000   0; 0  25000000]; % 4*4 matrix
S1 = [3.1683    3.1686    1.8716    1.8898]; % 1*4 vector
enforce_S_constraint=true;

% Create a struct so the parameters can be easily captured by the anonymous function
params=struct;
[params.N,params.K,params.M,params.C_l,params.H,params.A,params.C,...
    params.P,params.S1,params.enforce_S_constraint]=deal(...
    N,K,M,C_l,H,A,C,P,S1,enforce_S_constraint);
params_pass1=params;
params_pass1.enforce_S_constraint=false;

% First try without enforcing S2>=S1, then use the result as initial guess
B_initial_guess=zeros(N,M);
B_pass1=fminsearch(@(B) calculate_cost(B,params_pass1),B_initial_guess);
B_pass2=fminsearch(@(B) calculate_cost(B,params      ),B_pass1);

% This is only relevant once you find a valid solution
[~,S2]=calculate_cost(B_pass2,params);

function [cost,S2]=calculate_cost(B,params)
%calculate S2 from B and the other given inputs
[N,K,C_l,H,A,C,P,S1,enforce_S_constraint]=deal(...
    params.N,params.K,params.C_l,params.H,params.A,params.C,...
    params.P,params.S1,params.enforce_S_constraint);
d=zeros(N,1);
for j=1:1:N
    d(j) = (B(j,:)*P*B(j,:)')/((2^(2*C_l))-(norm(A(:,j))^2));
end
D_d = diag(d);

V_d=zeros(1,K);
sigma_d=zeros(1,K);
S2 = zeros(1,K); % intial value
for i=1:1:K   
    V_d(i)=C(i,:)*P*B'*H(i,:)'*inv(1+H(i,:)*(A'*D_d*A+B*P*B')*H(i,:)');
    sigma_d(i)=norm((V_d(i)*H(i,:)*B-C(i,:))*(P^(1/2)))^2+(V_d(i)^2)*(1+H(i,:)*A'*D_d*A*H(i,:)');
    
    S2(i)=0.5*log2((P(1,1))/sigma_d(:,i));   
end

if enforce_S_constraint && any(S1>S2)
    cost=inf;%set the cost to inf if the constraint is not satisfied
else
    cost=sum( (S2-S1).^2 );
end
end

[Adi Nor]

@Rik,Thanks for your help. your code helps me

Answer 1

[Matt J]

Well, the way you would approach it is using fmincon as below. Using your initial guess of zeros(N,M), fmincon is unable to find a solution so, assuming a feasible solution does exist, you will need to choose a less arbitrary initial guess. I urge you to get familiar with the fmincon documentation.

B_Optimal=doOptimization

Converged to an infeasible point.

fmincon stopped because the size of the current step is less than
the value of the step size tolerance but constraints are not
satisfied to within the value of the constraint tolerance.

B_Optimal = 2×2

46.2078 -19.7949 8.1964 -3.4964

●

function B_Optimal=doOptimization

% Given system parameters:
N = 2;
K = 4;
M=2;
C_l = 4;
H = [0.1185    0.2811; 0.3550    0.8224; 0.3260    0.9644; 0.5333    0.6083]; % 4*2 matrix
A = [-2 1; -1 1]; % 2*2 matrix
C = [7 -3; 7 -3; -2 1; -2 1]; % 4*2 matrix
P = [25000000   0; 0  25000000]; % 4*4 matrix
S1 = [3.1683    3.1686    1.8716    1.8898]; % 1*4 vector

Binitial = zeros(N,M); % intial value

% How can we  optimize the value of B matrix to achieve our goal? 

 B_Optimal=fmincon( @(B)norm( S2(B)-S1).^2,  Binitial,[],[],[],[],[],[],@nlcon);

 
    function [c,ceq]=nlcon(B)
        
        ceq=[];
        
        c=S1-S2(B);
        
    end


    function out=S2(B)
    %calculate S2 from B and the other given inputs
        for j=N:-1:1
            d(j) = (B(j,:)*P*B(j,:)')/((2^(2*C_l))-(norm(A(:,j))^2));
        end
        D_d = diag(d);
        
        for i=K:-1:1   
            V_d(i)=C(i,:)*P*B'*H(i,:)'*inv(1+H(i,:)*(A'*D_d*A+B*P*B')*H(i,:)');
            sigma_d(i)=norm((V_d(i)*H(i,:)*B-C(i,:))*(P^(1/2)))^2+(V_d(i)^2)*(1+H(i,:)*A'*D_d*A*H(i,:)');
            
            out(i)=0.5*log2((P(1,1))/sigma_d(:,i));   
        end
    end

end

[Stephan]

i agree - i tried with ga and also ga was not able to find a feasible solution that satisfied the constraint:

myOptProb

function myOptProb
% Given system parameters:
N = 2;
K = 4;
M = 2;
C_l = 4;
H = [0.1185    0.2811; 0.3550    0.8224; 0.3260    0.9644; 0.5333    0.6083]; % 4*2 matrix
A = [-2 1; -1 1]; % 2*2 matrix
C = [7 -3; 7 -3; -2 1; -2 1]; % 4*2 matrix
P = [25000000   0; 0  25000000]; % 4*4 matrix
S1 = [3.1683    3.1686    1.8716    1.8898]; % 1*4 vector
S2 = zeros(1,K); % intial value
B0 = zeros(N,M); % intial value

% How can we  optimize the value of B matrix to achieve our goal? 

opts = optimoptions('ga','MaxGenerations',20000,'PopulationSize',5000,'PlotFcn', @gaplotbestf);
[B, fval, exitflag, output] = ga(@fun,numel(B0),[],[],[],[],[],[],@(x)nonlcon(x,S1,S2),opts)

fun(B)

%calculate S2 from B and the other given inputs
function F = fun(B)
    B = reshape(B,N,M);
for j=1:1:N
    d(j) = (B(j,:)*P*B(j,:)')/((2^(2*C_l))-(norm(A(:,j))^2));
end
D_d = diag(d);

for i=1:1:K   
    V_d(i)=C(i,:)*P*B'*H(i,:)'*inv(1+H(i,:)*(A'*D_d*A+B*P*B')*H(i,:)');
    sigma_d(i)=norm((V_d(i)*H(i,:)*B-C(i,:))*(P^(1/2)))^2+(V_d(i)^2)*(1+H(i,:)*A'*D_d*A*H(i,:)');
    
    S2(i)=0.5*log2((P(1,1))/sigma_d(:,i));   
end

F = sum((S2-S1).^2);
end


function [c,ceq] = nonlcon(~,S1,S2)
c = S1(:)-S2(:);    % Compute nonlinear inequalities at x.
ceq = [];  % Compute nonlinear equalities at x.
end

end

ends up with:

Optimization terminated: no feasible point found.

B =

    0.5862    0.1652    1.4282   -1.0264


fval =

   54.7745


exitflag =

    -2


output = 

  struct with fields:

      problemtype: 'nonlinearconstr'
         rngstate: [1×1 struct]
      generations: 1
        funccount: 252512
          message: 'Optimization terminated: no feasible point found.'
    maxconstraint: 3.1686

[Rik]

How certain is it that such a B actually exists? It is not my field, so algebraically determining is such a matrix is possible in the first place is beyond my skills. However, I don't see anything to indicate that it would be impossible to find this out.