“Exiting: One or more of the residuals, duality ga

Hi All When I am running my code ,I got this Error message"

“Exiting: One or more of the residuals, duality gap, or total relative error has grown 100000 times greater than its minimum value so far: the primal appears to be infeasible (and the dual unbounded).(The dual residual < TolFun=1.00e-008.)”

Could you please help me how to solve this error

I look forward to hearing from You

Regards Bestun

[Geoff]

You could start by posting your code, or at least the optimisation function you are calling.

[Bestun]

% Example 1: dlo2(9, 6, 0, 0, 2, 5, 1) (Retaining wall; exact = (2 + pi) / 2) %

% %

function dlo2(xmax, ymax, edgeA, edgeB, edgeC, edgeD, phiDegs);

nodes = createNodes(xmax, ymax);

discs = createDiscontinuities(nodes, xmax, ymax, edgeA, edgeB, edgeC,edgeD);

B = compatibilityMatrix(nodes, discs);

[objP padN N] = plasticMultiplierTerms(nodes, discs, phiDegs);

fD = selfWeight(nodes, discs, ymax);

fL = unitLoad(discs);

tied = tieDiscs(discs);

[vars, soln] = solve(B, N, padN, fD, fL, tied, objP);

plotMechanism(nodes, discs, vars, xmax, ymax);

soln

%===================================================================================%

% SETUP NODES %

function nodes = createNodes(nx, ny)

nodes = [];

for y = 0 : ny

for x = 0 : nx

nodes = [nodes; size(nodes, 1) + 1 x y];

end

end

% SETUP DISCONTINUITY ELEMENTS %

function discs = createDiscontinuities(nodes, xmax, ymax, edgeA, edgeB, edgeC, edgeD)

discs = [];

edgeA = sparse(3, 1) + edgeA; edgeB = sparse(3, 1) + edgeB; edgeC = sparse(3, 1) + edgeC; edgeD = sparse(3, 1) + edgeD;

for i = 1 : size(nodes, 1) - 1

for j = i + 1 : size(nodes, 1)

edgeType = edgeA(1) * ((nodes(i, 3) + nodes(j, 3) == 0) && (nodes(i, 2) <= (xmax - edgeA(3))) && (nodes(j, 2) <= (xmax - edgeA(3))))...

+ edgeA(2) * ((nodes(i, 3) + nodes(j, 3) == 0) && (nodes(i, 2) >= (xmax - edgeA(3))) && (nodes(j, 2) >= (xmax - edgeA(3))))...

+ edgeB(1) * (nodes(i, 2) + nodes(j, 2) == 2 * xmax) ...

+ edgeC(1) * ((nodes(i, 3) + nodes(j, 3) == 2 * ymax) && (nodes(i, 2) <= (xmax - edgeC(3))) && (nodes(j, 2) <= (xmax - edgeC(3))))...

+ edgeC(2) * ((nodes(i, 3) + nodes(j, 3) == 2 * ymax) && (nodes(i, 2) >= (xmax - edgeC(3))) && (nodes(j, 2) >= (xmax - edgeC(3))))...

+ edgeD(1) * (nodes(i, 2) + nodes(j, 2) == 0);

if (gcd(abs(nodes(j, 2) - nodes(i, 2)), abs(nodes(j, 3) - nodes(i, 3))) < 1.0001)

discs = [discs; i j edgeType sqrt((nodes(i, 2) - nodes(j, 2))^2 + (nodes(i, 3) - nodes(j, 3))^2)];

end

end

end

% SETUP COMPATIBILIY MATRIX %

function B = compatibilityMatrix(nodes, discs)

[nCons nVars] = deal(2 * size(nodes, 1), 2 * size(discs, 1));

B = sparse(nCons, nVars);

for i = 1 : size(discs, 1)

[n1 n2 len] = deal(discs(i, 1), discs(i, 2), discs(i, 4));

[conN1 conN2 var] = deal(2 * nodes(n1, 1) - 1, 2 * nodes(n2, 1) - 1, 2 * i - 1);

[cosine sine] = deal((nodes(n1, 2) - nodes(n2, 2)) / len, (nodes(n1, 3) - nodes(n2, 3)) / len);

[cosine sine];

B_local = [cosine -sine; sine cosine];

B([conN1 conN1 + 1], [var var + 1]) = B_local;

B([conN2 conN2 + 1], [var var + 1]) = -B_local;

end

% SETUP PLASTIC MULTIPLIER MATRIX %

function [objP, padN, N] = plasticMultiplierTerms(nodes, discs, phiDegrees)

objP = sparse(0, 0); padN = sparse(0, 0); N = sparse(0, 0); count = 1; tan_phi = tan(pi * phiDegrees / 180); Ce = 0; %z =1;

for i = 1 : size(discs, 1)

if (discs(i, 3) < 2) % i.e. if rigid or symmetry

[eff_Ce eff_tanPhi] = deal(Ce * (discs(i, 3) ~= 1), tan_phi * (discs(i, 3) ~= 1));

[eff_Ce eff_tanPhi];

[con var] = deal(2 * count - 1, 2 * count - 1);

N_local = [1 -1; eff_tanPhi eff_tanPhi];

N([con con + 1], [var var + 1]) = N_local;

padN_local = sparse(2, 2 * size(discs, 1));

padN_local(1, 2 * i - 1) = -1;

padN_local(2, 2 * i) = -1;

padN = [padN; padN_local];

[n1 n2 len] = deal(discs(i, 1), discs(i, 2), discs(i, 4));

x1 = nodes(n1,2);

y1 = nodes(n1,3);

x2 = nodes(n2,2);

y2 = nodes(n2,3);

theta = atan((y2-y1)/(x2-x1));

theta = 180 *theta/pi;

[sine] = deal(nodes(n1, 3) - nodes(n2, 3)) / len; %Why value of sine is negative for Discs 4

[rw a Hw Yw rd]= deal(9.81, 0.017, 0, -1000, 8.8);

if theta == 0 && y1 >= Yw && y2 >= Yw; % above water table

Ce = rw *(y1-Yw)*exp(a*rw*Yw)*exp(-a*rw*y1)*tan_phi*(x2-x1); %OKKKKKK % this case is already neglected since we used this condition... if (discs(i, 3) < 2): this case acts on discs 6 (for 1*1 square) , on this discs apparent cohesion is neglected

elseif theta == 0 && y1 < Yw && y2 < Yw; % under water table

Ce = rw *(y1-Yw)*tan_phi*(x2-x1); %OKKKKKK

end

if theta == 90 && y1 >= Yw && y2 > Yw ; % above water table

Ce = ((-1)*tan_phi/((a)^2*rw))*((1-a*rw*Yw + a*rw*y2)*exp(a*rw*(Yw-y2))-(1-a*rw*Yw + a*rw*y1)*exp(a*rw*(Yw-y1))); %OKKKKKK

elseif theta == 90 && y1 < Yw && y2 <= Yw; % under water table

Ce = rw*tan_phi*(((y2^2/2)-Yw*y2)-((y1^2/2)-Yw*y1)); %OKKKKKK

elseif theta == 90 && y1 <= Yw && y2>Yw; % case is between

Ce = rw*tan_phi* (((Yw^2/2)-Yw^2)-((y1^2/2)-Yw*y1))-(tan_phi/((a)^2*rw))*((1-a*rw*Yw + a*rw*y2)*exp(a*rw*(Yw-y2))-1); %OKKKKKK

end

if theta ~= 0 && theta ~= 90 && y1 >= Yw && y2 > Yw ; % above water table;

Ce =(((-1)*exp(a*rw*(Yw-y1))*tan_phi)/((a^2)*rw*sine))*(exp((-1)*a*rw*discs(i,4)*sine)*(a*rw*(discs(i,4)*sine-(Yw-y1))+1)+ (a*rw*(Yw-y1))-1); %OKKKKKK

elseif theta ~= 0 && theta ~= 90 && y1 < Yw && y2 <= Yw; % under water table

Ce = rw*tan_phi*(((discs(i,4)^2/2)*sine)+ discs(i,4)*(y1-Yw)); %OKKKKKK % Why the program takes the biggest value of Apparent cohesion

elseif theta ~= 0 && theta ~= 90 && y1 < Yw && y2>Yw; L1 = 0; L2=0 ; % case is between

L1= abs((Yw-y1)* sine); % here I used abs values for L1 since all values of sine are negative

L2 = abs((y2-Yw)*sine); % here I used abs values for L1 since all values of sine are negative

Ce = rw*tan_phi*(((L1^2/2)*sine)+ L1*(y1-Yw))-((tan_phi)/((a^2)*rw*sine))*(((1+ a*rw*L2*sine)*exp(-a*rw*L2*sine))-((1+ a*rw*Yw*sine)*exp(-a*rw*Yw*sine))); %OKKKKKK

% z

end

objP = [objP Ce Ce];

count = count + 1;

end

% z = z+1;

end

% SELF WEIGHT LOADS %

function fD = selfWeight(nodes, discs, ymax)

fD = []; x = 1; stripWeight = 0;

for i = 1 : size(discs, 1)

[n1 n2 len] = deal(discs(i, 1), discs(i, 2), discs(i, 4));

x1 = nodes(n1,2);

y1 = nodes(n1,3);

x2 = nodes(n2,2);

y2 = nodes(n2,3);

if x1 ~= x2

theta = atan((y2-y1)/(x2-x1));

theta = 180 *theta/pi;

[cosine sine] = deal((nodes(n1, 2) - nodes(n2, 2)) / len, (nodes(n1, 3) - nodes(n2, 3)) / len);

[cosine sine];

rd = 0; rsat =0.0; a = 0; Hw =0; Yw =0; rw =0;

[rw a Hw Yw rd rsat]= deal(9.81, 0.017, 0, -1000, 8.8, 15.2);

if theta ~= 0 && theta ~= 90 && y1 >= Yw && y2 > Yw ; % above water table;

stripWeight =(cosine/(a*rw))*(((a*rw*rd*ymax)-((rsat-rd)*exp(a*rw*Yw)*exp((-1)*a*rw*ymax))-(a*rw*rd*y1))*discs(i,4)-(a*rw*rd*sine*(discs(i,4))^2/2)-((rsat-rd)*exp(a*rw*(Yw-y1))*exp((-1)*a*rw*discs(i,4)*sine)/(sine*a*rw)) +(((rsat-rd)*exp(a*rw*(Yw-y1)))/(sine*a*rw))); %OKKKKK

x;

elseif theta ~= 0 && theta ~= 90 && y1 < Yw && y2 <= Yw && Yw < ymax; % under water table

stripWeight = ((rsat*Yw*discs(i,4))-((discs(i,4)^2/2)*rsat*sine)-(rsat*y1*discs(i,4)))*cosine...

+ ((rd *ymax -((rsat-rd)/a*rw)*exp(a*rw*Yw)*exp((-1)*rw*ymax)+(rsat-rd)/a*rw))*(x2-x1); % not checked In this case "WE NEED 1*2 RECTANGULAR GEOMETRY"

elseif theta ~= 0 && theta ~= 90 && y1 < Yw && y2 <= Yw && Yw == ymax; % under water table

stripWeight = ((rsat*Yw*discs(i,4))-((discs(i,4)^2/2)*rsat*sine)-(rsat*y1*discs(i,4)))*cosine; % OKKKKKKK

elseif theta ~= 0 && theta ~= 90 && y1 < Yw && y2>Yw; L1 = 0; L2=0 ;% case is between

L1= abs((Yw-y1)* sine); % here I used abs values for L1 since all values of SINE are negative

L2 = abs((y2-Yw)*sine); % here I used abs values for L1 since all values of SINE are negative

stripWeight = (cosine/(a*rw))*(((a*rw*rd*ymax)-((rsat-rd)*exp(a*rw*Yw)*exp((-1)*a*rw*ymax))-(a*rw*rd*y1))*L2-(a*rw*rd*sine*(L2)^2/2)-((rsat-rd)*exp((-1)*a*rw*L2*sine)/(sine*a*rw))+((rsat-rd)/(sine*a*rw)))...

+(rsat*Yw*L1-(L1^2/2)*rsat*sine-rsat*y1*L1)*cosine...

+(1/(a*rw))*(((a*rw*rd*ymax)-((rsat-rd)*exp(a*rw*Yw)*exp((-1)*a*rw*ymax))-(a*rw*rd*Yw))+ (rsat-rd))*(L1*cosine); %OKKKKKK

end

if theta == 0 && y1 <= Yw && y2 <= Yw && Yw ==ymax; % under water table

stripWeight = rsat*(ymax-y1)*(x2-x1); %OKKKKKK

elseif theta == 0 && y1 <= Yw && y2 <= Yw && Yw < ymax; % CASE IN BETWEEN

stripWeight = rsat*(Yw-y1)*(x2-x1)+(rd*ymax-((rsat-rd)*exp(a*rw*(Yw-ymax))/(a*rw))-rd*Yw +((rsat-rd)/(a*rw)))*(x2-x1); %OKKKKKKK

[stripWeight]

elseif theta == 0 && y1 > Yw && y2 > Yw ; % above water table

x1 = nodes(n1,2);

stripWeight = (rd*ymax -((rsat-rd)*exp(a*rw*(Yw-ymax))/(a*rw))-rd*y1 + ((rsat-rd)*exp(a*rw*(Yw-y1))/(a*rw)))*(x2-x1); %OKKKKKK

end

end

fD = [fD; -sine * stripWeight; -cosine * stripWeight];

x =x+1;

end

% EXTERNAL LIVE LOADS %

function fL = unitLoad(discs);

fL = sparse(1, 2 * size(discs, 1));

for i = 1 : size(discs, 1)

if (discs(i, 3) >= 3)

fL(1, 2 * i) = 1; %normal stress

end

end

% LINK ELEMENT DISPLACEMENTS (ONLY IF RIGID LOADS) %

function tied = tieDiscs(discs); % for simplicity implement by adding constraints

tied = sparse(0, 2 * size(discs, 1)); ilast = 0;

for i = 1 : size(discs, 1)

if (discs(i, 3) >= 4) % i.e. if rigid load

if ((ilast ~= 0) || (discs(i, 3) == 5))

extra = sparse(2, 2 * size(discs, 1));

extra(1, 2 * i - 1) = 1;

if ilast ~= 0

extra(1, 2 * ilast - 1) = -1;

extra(2, 2 * i) = 1; extra(2, 2 * ilast) = -1;

end

tied = [tied; extra];

end

ilast = i;

end

end

% FORM AND SOLVE LINEAR PROGRAMMING PROBLEM %

function [variables, solution] = solve(B, N, padN, fD, fL, tied, objP);

bigNumber = 1000; warning off MATLAB:divideByZero;

[padB padFL padTied] = deal(sparse(size(B, 1), size(N, 2)), sparse(1, size(N, 2)), sparse(size(tied, 1), size(N, 2)));

equalityMatrix = [B padB; padN N; fL padFL; tied padTied];

equalityRHS = sparse(size(equalityMatrix, 1), 1);

equalityRHS(size(B, 1) + size(N, 1) + 1, 1) = 1;

obj = [fD; objP']; % objP' contains plastic multiplier terms

lowerB = [-bigNumber * ones(size(B, 2), 1); sparse(size(N, 2), 1)];

upperB = [bigNumber * ones(size(B, 2), 1); bigNumber * ones(size(N, 2), 1)];

options=optimset('Display','iter','Tolfun',0.1)

[variables, solution] = linprog(obj, [], [], equalityMatrix, equalityRHS, lowerB, upperB);

% PLOT MECHANISM %

function plotMechanism(nodes, discs, variables, xmax, ymax);

clf; hold on; axis equal;

set(gca,'XLim',[0 xmax], 'XTick', [0:1:xmax], 'YLim',[0 ymax], 'YTick', [0:1:ymax], 'FontSize', 12)

xlabel('x', 'FontSize', 14); ylabel('y', 'FontSize', 12);

for pass = 1 : 2 % plot active elements on top of inactive elements

for i = 1 : size(discs, 1)

[xp yp] = deal([nodes(discs(i, 1), 2); nodes(discs(i, 2), 2)], [nodes(discs(i, 1), 3); nodes(discs(i, 2), 3)]);

if ((pass == 2) && ((abs(variables(2 * i)) > 0.001) || (abs(variables((2 * i) - 1)) > 0.001)) && (discs(i, 3) ~= 1))

plot(xp, yp, 'r-', 'LineWidth', 3); % active = red

elseif (pass == 1)

plot(xp, yp, '-', 'Color', [0.9 0.9 0.9]); % inactive = light grey

end

end

end

Answer 1

[Bestun]

        Thanks Geoff for reply
The above is my code.
When I tried to use these parameters 
dlo2(1, 1, 0, 0, 2, 5, 36)
I got the above error ...could you please tell me what is the wrong

regards Bestun