Good afternoon,

Introduction

I need to solve the following optimization problem, whose theoretical description can be found in sections 7.7.1 and 7.7.2 of this link.

Definition of fmincon and the target function

function [v] = LP_Irr_LDPC(N,Ebn0)

options = optimoptions('fmincon','Display','iter','Algorithm','interior-point','MaxIter', 10000, 'MaxFunEvals', 10000, 'StepTolerance', 1e-100);
fun = @(v) -sum(v(1:N)./(1:N));

A = [];
b = [];
Aeq = [0, ones(1,N-1)];
beq = 1;
lb = zeros(1,N);
ub = [0, ones(1,N-1)];
nonlcon = @(v)DensEv_BEC(v,Ebn0);
l0 = [0 rand(1,N-1)];
l0 = round(l0./sum(l0),2);

v = fmincon(fun,l0,A,b,Aeq,beq,lb,ub,nonlcon,options);
end

Restrictions

function [c, ceq] = DensEv_BEC(v,Ebn0)

% It is also needed to modify this function, as you cannot pass parameters
% from others to it. 

h = [0 rand(1,999)];
h = h./sum(h);

syms x;
X = x.^(0:(length(h)-1));
R = h*transpose(X);

ebn0 = 10^(Ebn0/10);
Rm = 1;
LLR = (-50:50);
p03 = 0.3;
LLR03 = log((1-p03)/p03);
r03 = 1 - p03;
noise03 = (2*r03*Rm*ebn0)^-1;
pf03 = normpdf(LLR, LLR03, noise03);
sumpf03 = sum(pf03(1:length(pf03)/2));

divisions = 200;

Aj = zeros(1, divisions);
rho = zeros(1, divisions);
xj = zeros(1, divisions);
N = 2000; % Length(v) -> Same value as in 'Complete.m'

for j=1:1:divisions
    xj(j) = sumpf03*j/divisions;
    rho(j) = subs(R,x,1-xj(j));
    Aj(j) = 1 - rho(j);
end

c = zeros(1, length(xj));
lambda = zeros(1, length(Aj));
for j = 1:1:length(xj)
    lambda(j) = sum(v(2:N).*(Aj(j).^(1:(N-1))));
    c(j) = sumpf03*lambda(j) - xj(j);
end

save Almacen
%ceq = [];
ceq = sum(v)-1;
end

Main function

The result of the above optimization problem has to be used within the next script. However, as the optimization problem fails to converge to any value that satisfies all conditions (namely, that the sum of elements of v is equal to the unit) the MacKayNeal function.m gives error. You can find more information on how this function works in the following link within this same page, where I ask another unresolved question about it:

clear;clc

N = 2000;
r = 0.5;
%EbN0 = [0 2 4 6 8 10 12]; % EbN0 in dB
EbN0 = 1;
iterations = 15; % Número de iteraciones del codificador (suave)
frames = 5; % Number of frames (N bits per frame)

for i = 1:length(EbN0)
    % LP - Distribuciones de bits por filas y por columnas
    v = LP_Irr_LDPC(N,EbN0(i));
    load('Almacen.mat');
    % DE - Distribuciones de bits por filas y por columnas
    %[v_DE h_DE]

    % LP - Determinación de las matrices generadora y de paridad
    H_LP = MacKayNeal(N,r,v,h);
    G_LP = GaussElm(H_LP,N,N*r);
    % DE - Determinación de las matrices generadora y de paridad
    %H_DE = MacKayNeal(N,r,v_DE,h_DE);
    %G_DE = GaussElm(H_DE,N,N*r);

    BER_LP(i) = 0;
    %BER_DE(i) = 0;
    for j = 1:frames
        % Encoding message
        u = transpose(round(rand(N*r,1)));       
        tx_LP = u*G_LP;
        %TX_DE = u*G_DE;

        % BPSK modulation
        bpskMod_LP = 2*tx_LP - 1;
        %bpskMod_DE = 2*tx_DE - 1;

        % Addition of AWGN
        N0 = 1/(10^(EbN0(i)/10));
        rx_LP = bpskMod_LP + sqrt(N0/2)*randn(size(bpskMod_LP));
        %rx_DE = bpskMod_DE + sqrt(N0/2)*randn(size(bpskMod_DE));

        % Decoding
        vhat_LP = decodeLogDomain(rx_LP, H_LP, N0, iterations); % Hay que hacer ésto
        %vhat_DE = decodeLogDomain(rx_DE, H_DE, N0, iterations);

        % Get BER
        [num_LP, rat_LP] = biterr(vhat_LP, tx_LP);
        BER_LP(i) = BER_LP(i)+rat_LP;
        %[num_DE, rat_DE] = biterr(vhat_DE, tx_DE);
        %BER_DE(i) = BER_DE(i)+rat_DE;
    end

    % Get average of BER
    BER_LP(i) = BER_LP(i)/frames;
    %BER_DE(i) = BER_DE(i)/frames;
end

% Plot the result
semilogy(EbN0, BER_LP, 'o-');
grid on; hold on;
%semilogy(EbN0, BER_DE, 'x-');
%grid on; hold off;