MATLAB: Nonlinear inequality constraint in Fmincon produces optimized x values only where constraint value equals zero

fminconnonlinear inequality constraint

When using the nonlinear inequality constraint function in fmincon, the optimized x values have been limited where c(x) = 0 instead of c(x) <= 0.
Objective function
 function fun = fgoptprac(x)
%Load these initial values into the workspace before running (x0, lb, ub)
%x0 = [.75;.05];
%lb = [.6;.02];
%ub = [.9;.1];
FGcost = 3;
PCcost = 120;
FAcost = 50;
FGion = 0.6245;
PCion = 0.0466;
FAion = 0.0408;
Md = 88; 
Vs = 62.4; 
As = 89.7; 
C0 = 2.52;
a = [1.5323;3.2456;4.4511;5.4375;1.3112;5.9278;3.5471;9.7522;1.1443;10.1803;11.6411;13.5043;10.6609;10.8304;9.7646;3.0424;8.7430];
b = [-36.4593;-525.0273;-69.7738;-271.1375;-108.3028;-89.6767;51.8115;-196.8059;-1631.6804;476.9985;-627.1805;-60.4963;-482.9669;-804.7777;-1356.6831;-1184.5658;-2072.8314];
c = [7.6043;2.3144;-20.9358;-22.1573;-53.1397;-20.613;-38.6912;-9.8625;-32.7446;-15.3286;4.9522;20.7256;1.1675;8.0776;10.1853;-12.5374;-2.0985];
d = [42.5919;578.303;76.8362;311.1605;153.3032;77.0661;-58.4724;189.645;1864.9115;-566.1635;675.7103;30.411;509.1583;881.4849;1492.1294;1378.4112;2299.2668];
e = [-9.8215;-9.7098;38.1797;35.9569;93.4644;25.7382;61.4931;4.5864;49.8227;24.3662;-17.3297;-44.1011;-12.3811;-28.7089;-36.184;33.7396;-18.6562];
f = [40.2246;607.0888;105.9604;285.1265;93.4722;195.0361;-19.3576;263.1551;1827.219;-543.1824;691.8505;65.2143;582.846;906.2586;1555.2136;1177.1889;2392.7099];
texp = [.08;.29;1;2;3;4;5;7;10;14;21;28;35;42;49;63;77];
t0 = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1]; 
    function C = conc(fgoptprac)
        C = (((x(1)*FGion +x(2)*PCion + (1 - x(1) - x(2))*FAion)*Md*1000/Vs) - C0);
    end
    function Cb = Bulk(fgoptprac)
        Cb = ((x(1)*FGion +x(2)*PCion + (1 - x(1) - x(2))*FAion)*Md*1000/Vs);
    end
    function cost = COST(fgoptprac)
        cost = (x(1)*FGcost + x(2)*PCcost + (1 - x(1) - x(2))*FAcost + 5);
    end
    function ionr = IONR(fgoptprac)
        ionr = (x(1)*FGion +x(2)*PCion + (1 - x(1) - x(2))*FAion);
    end
    function delta = DELTA(fgoptprac)
        delta = 1 - (Mtc*As/(Bulk*Vs));
        function mtc = Mtc(DELTA)
         k1 = (conc*t0);
         k2 = (conc*sqrt(texp));
         k3 = (conc*texp);
         X = [k1 k2 k3];
         Mint = [a(1)*x(1) + b(1)*x(2) + c(1)*(1 - x(1) - x(2)) + d(1)*x(1)*x(2) + ... 
         e(1)*x(1)*(1-x(1)-x(2)) + f(1)*x(2)*(1-x(1)-x(2)); 
         a(2)*x(1) + b(2)*x(2) + c(2)*(1 - x(1) - x(2)) + d(2)*x(1)*x(2) + ...
         e(2)*x(1)*(1-x(1)-x(2)) + f(2)*x(2)*(1-x(1)-x(2));
         a(3)*x(1) + b(3)*x(2) + c(3)*(1 - x(1) - x(2)) + d(3)*x(1)*x(2) + ...
         e(3)*x(1)*(1-x(1)-x(2)) + f(3)*x(2)*(1-x(1)-x(2));
         a(4)*x(1) + b(4)*x(2) + c(4)*(1 - x(1) - x(2)) + d(4)*x(1)*x(2) + ...
         e(4)*x(1)*(1-x(1)-x(2)) + f(4)*x(2)*(1-x(1)-x(2));
         a(5)*x(1) + b(5)*x(2) + c(5)*(1 - x(1) - x(2)) + d(5)*x(1)*x(2) + ...
         e(5)*x(1)*(1-x(1)-x(2)) + f(5)*x(2)*(1-x(1)-x(2));
         a(6)*x(1) + b(6)*x(2) + c(6)*(1 - x(1) - x(2)) + d(6)*x(1)*x(2) + ...
         e(6)*x(1)*(1-x(1)-x(2)) + f(6)*x(2)*(1-x(1)-x(2));
         a(7)*x(1) + b(7)*x(2) + c(7)*(1 - x(1) - x(2)) + d(7)*x(1)*x(2) + ...
         e(7)*x(1)*(1-x(1)-x(2)) + f(7)*x(2)*(1-x(1)-x(2));
         a(8)*x(1) + b(8)*x(2) + c(8)*(1 - x(1) - x(2)) + d(8)*x(1)*x(2) + ...
         e(8)*x(1)*(1-x(1)-x(2)) + f(8)*x(2)*(1-x(1)-x(2));
         a(9)*x(1) + b(9)*x(2) + c(9)*(1 - x(1) - x(2)) + d(9)*x(1)*x(2) + ...
         e(9)*x(1)*(1-x(1)-x(2)) + f(9)*x(2)*(1-x(1)-x(2));
         a(10)*x(1) + b(10)*x(2) + c(10)*(1 - x(1) - x(2)) + d(10)*x(1)*x(2) + ...
         e(10)*x(1)*(1-x(1)-x(2)) + f(10)*x(2)*(1-x(1)-x(2));
         a(11)*x(1) + b(11)*x(2) + c(11)*(1 - x(1) - x(2)) + d(11)*x(1)*x(2) + ...
         e(11)*x(1)*(1-x(1)-x(2)) + f(11)*x(2)*(1-x(1)-x(2));
         a(12)*x(1) + b(12)*x(2) + c(12)*(1 - x(1) - x(2)) + d(12)*x(1)*x(2) + ...
         e(12)*x(1)*(1-x(1)-x(2)) + f(12)*x(2)*(1-x(1)-x(2));
         a(13)*x(1) + b(13)*x(2) + c(13)*(1 - x(1) - x(2)) + d(13)*x(1)*x(2) + ...
         e(13)*x(1)*(1-x(1)-x(2)) + f(13)*x(2)*(1-x(1)-x(2));
         a(14)*x(1) + b(14)*x(2) + c(14)*(1 - x(1) - x(2)) + d(14)*x(1)*x(2) + ...
         e(14)*x(1)*(1-x(1)-x(2)) + f(14)*x(2)*(1-x(1)-x(2));
         a(15)*x(1) + b(15)*x(2) + c(15)*(1 - x(1) - x(2)) + d(15)*x(1)*x(2) + ...
         e(15)*x(1)*(1-x(1)-x(2)) + f(15)*x(2)*(1-x(1)-x(2));
         a(16)*x(1) + b(16)*x(2) + c(16)*(1 - x(1) - x(2)) + d(16)*x(1)*x(2) + ...
         e(16)*x(1)*(1-x(1)-x(2)) + f(16)*x(2)*(1-x(1)-x(2));
         a(17)*x(1) + b(17)*x(2) + c(17)*(1 - x(1) - x(2)) + d(17)*x(1)*x(2) + ...
         e(17)*x(1)*(1-x(1)-x(2)) + f(17)*x(2)*(1-x(1)-x(2));];
         y = cumsum(Mint);
         beta = regress(y,X);
         K1 = beta(1);
         K2 = beta(2);
         K3 = beta(3);
         tc = (K2/(2*K3))^2;
         mtc = (K1*conc + K2*conc*sqrt(tc) + K3*conc*tc);
         end
    end
fun = (COST/(DELTA*IONR));
%Insert following code [x,fval]=.... directly into command window
%[x,fval] = fmincon(@fgoptprac,x0,[],[],[],[],lb,ub,@confun);
end
Constraint function
    function [con,ceq]=confun(x)
FGcost = 3;
PCcost = 120;
FAcost = 50;
FGion = 0.6245; 
PCion = 0.0466;
FAion = 0.0408;
Md = 88;
Vs = 62.4;
As = 89.7;
C0 = 2.52;
a = [1.5323;3.2456;4.4511;5.4375;1.3112;5.9278;3.5471;9.7522;1.1443;10.1803;11.6411;13.5043;10.6609;10.8304;9.7646;3.0424;8.7430];
b = [-36.4593;-525.0273;-69.7738;-271.1375;-108.3028;-89.6767;51.8115;-196.8059;-1631.6804;476.9985;-627.1805;-60.4963;-482.9669;-804.7777;-1356.6831;-1184.5658;-2072.8314];
c = [7.6043;2.3144;-20.9358;-22.1573;-53.1397;-20.613;-38.6912;-9.8625;-32.7446;-15.3286;4.9522;20.7256;1.1675;8.0776;10.1853;-12.5374;-2.0985];
d = [42.5919;578.303;76.8362;311.1605;153.3032;77.0661;-58.4724;189.645;1864.9115;-566.1635;675.7103;30.411;509.1583;881.4849;1492.1294;1378.4112;2299.2668];
e = [-9.8215;-9.7098;38.1797;35.9569;93.4644;25.7382;61.4931;4.5864;49.8227;24.3662;-17.3297;-44.1011;-12.3811;-28.7089;-36.184;33.7396;-18.6562];
f = [40.2246;607.0888;105.9604;285.1265;93.4722;195.0361;-19.3576;263.1551;1827.219;-543.1824;691.8505;65.2143;582.846;906.2586;1555.2136;1177.1889;2392.7099];
texp = [.08;.29;1;2;3;4;5;7;10;14;21;28;35;42;49;63;77];
t0 = [1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1];
    function C = conc2(confun)
        C = (((x(1)*FGion +x(2)*PCion + (1 - x(1) - x(2))*FAion)*Md*1000/Vs) - C0);
    end
    function Cb = Bulk2(confun)
        Cb = ((x(1)*FGion +x(2)*PCion + (1 - x(1) - x(2))*FAion)*Md*1000/Vs);
    end
    function cost = COST2(confun)
        cost = (x(1)*FGcost + x(2)*PCcost + (1 - x(1) - x(2))*FAcost + 5);
    end
    function ionr = IONR2(confun)
        ionr = (x(1)*FGion +x(2)*PCion + (1 - x(1) - x(2))*FAion);
    end
    function delta = DELTA2(confun)
        delta = 1 - (Mtc*As/(Bulk2*Vs));
        function mtc = Mtc(DELTA2)
         k1 = (conc2*t0);
         k2 = (conc2*sqrt(texp));
         k3 = (conc2*texp);
         X = [k1 k2 k3];
         Mint = [a(1)*x(1) + b(1)*x(2) + c(1)*(1 - x(1) - x(2)) + d(1)*x(1)*x(2) + ... 
         e(1)*x(1)*(1-x(1)-x(2)) + f(1)*x(2)*(1-x(1)-x(2)); 
         a(2)*x(1) + b(2)*x(2) + c(2)*(1 - x(1) - x(2)) + d(2)*x(1)*x(2) + ...
         e(2)*x(1)*(1-x(1)-x(2)) + f(2)*x(2)*(1-x(1)-x(2));
         a(3)*x(1) + b(3)*x(2) + c(3)*(1 - x(1) - x(2)) + d(3)*x(1)*x(2) + ...
         e(3)*x(1)*(1-x(1)-x(2)) + f(3)*x(2)*(1-x(1)-x(2));
         a(4)*x(1) + b(4)*x(2) + c(4)*(1 - x(1) - x(2)) + d(4)*x(1)*x(2) + ...
         e(4)*x(1)*(1-x(1)-x(2)) + f(4)*x(2)*(1-x(1)-x(2));
         a(5)*x(1) + b(5)*x(2) + c(5)*(1 - x(1) - x(2)) + d(5)*x(1)*x(2) + ...
         e(5)*x(1)*(1-x(1)-x(2)) + f(5)*x(2)*(1-x(1)-x(2));
         a(6)*x(1) + b(6)*x(2) + c(6)*(1 - x(1) - x(2)) + d(6)*x(1)*x(2) + ...
         e(6)*x(1)*(1-x(1)-x(2)) + f(6)*x(2)*(1-x(1)-x(2));
         a(7)*x(1) + b(7)*x(2) + c(7)*(1 - x(1) - x(2)) + d(7)*x(1)*x(2) + ...
         e(7)*x(1)*(1-x(1)-x(2)) + f(7)*x(2)*(1-x(1)-x(2));
         a(8)*x(1) + b(8)*x(2) + c(8)*(1 - x(1) - x(2)) + d(8)*x(1)*x(2) + ...
         e(8)*x(1)*(1-x(1)-x(2)) + f(8)*x(2)*(1-x(1)-x(2));
         a(9)*x(1) + b(9)*x(2) + c(9)*(1 - x(1) - x(2)) + d(9)*x(1)*x(2) + ...
         e(9)*x(1)*(1-x(1)-x(2)) + f(9)*x(2)*(1-x(1)-x(2));
         a(10)*x(1) + b(10)*x(2) + c(10)*(1 - x(1) - x(2)) + d(10)*x(1)*x(2) + ...
         e(10)*x(1)*(1-x(1)-x(2)) + f(10)*x(2)*(1-x(1)-x(2));
         a(11)*x(1) + b(11)*x(2) + c(11)*(1 - x(1) - x(2)) + d(11)*x(1)*x(2) + ...
         e(11)*x(1)*(1-x(1)-x(2)) + f(11)*x(2)*(1-x(1)-x(2));
         a(12)*x(1) + b(12)*x(2) + c(12)*(1 - x(1) - x(2)) + d(12)*x(1)*x(2) + ...
         e(12)*x(1)*(1-x(1)-x(2)) + f(12)*x(2)*(1-x(1)-x(2));
         a(13)*x(1) + b(13)*x(2) + c(13)*(1 - x(1) - x(2)) + d(13)*x(1)*x(2) + ...
         e(13)*x(1)*(1-x(1)-x(2)) + f(13)*x(2)*(1-x(1)-x(2));
         a(14)*x(1) + b(14)*x(2) + c(14)*(1 - x(1) - x(2)) + d(14)*x(1)*x(2) + ...
         e(14)*x(1)*(1-x(1)-x(2)) + f(14)*x(2)*(1-x(1)-x(2));
         a(15)*x(1) + b(15)*x(2) + c(15)*(1 - x(1) - x(2)) + d(15)*x(1)*x(2) + ...
         e(15)*x(1)*(1-x(1)-x(2)) + f(15)*x(2)*(1-x(1)-x(2));
         a(16)*x(1) + b(16)*x(2) + c(16)*(1 - x(1) - x(2)) + d(16)*x(1)*x(2) + ...
         e(16)*x(1)*(1-x(1)-x(2)) + f(16)*x(2)*(1-x(1)-x(2));
         a(17)*x(1) + b(17)*x(2) + c(17)*(1 - x(1) - x(2)) + d(17)*x(1)*x(2) + ...
         e(17)*x(1)*(1-x(1)-x(2)) + f(17)*x(2)*(1-x(1)-x(2));];
         y = cumsum(Mint);
         beta = regress(y,X);
         K1 = beta(1);
         K2 = beta(2);
         K3 = beta(3);
         tc = (K2/(2*K3))^2;
         mtc = (K1*conc2 + K2*conc2*sqrt(tc) + K3*conc2*tc);
         end
    end
        a00 = -208.28;
        a01 = 815.50;
        a02 = -1009.08;
        a03 = 404.27;
        a10 = 3082.11;
        a11 = -11411.57;
        a12 = 13930.50;
        a13 = -5622.92;
        a20 = -2418.38;
        a21 = 8758.37;
        a22 = -10446.88;
        a23 = 4116.15;
    function Sc = Cstr(confun)
        Sc = a00 + a01*x(1) + a02*(x(1))^2 + a03*(x(1))^3 + a10*(x(2)*10) ...
            + a11*(x(2)*10)*x(1) + a12*(x(2)*10)*(x(1))^2 + a13*(x(2)*10)*(x(1))^3 ...
            + a20*(x(2)*10)^2 + a21*(x(2)*10)^2*x(1) ...
            + a22*(x(2)*10)^2*(x(1))^2 + a23*(x(2)*10)^2*(x(1))^3;
    end
con(1) = -DELTA2;
con(2) = -Cstr + 4.137;
ceq = [];
end
The output of the code above has been confirmed to be correct when run with a defined x. However, when the fmincon procedure is called with these constraint functions, the output x is given when con(2) = -Cstr + 4.137 = 0.
Effectively this only outputs values of x where con(2) = 0 instead of minimizing the function, fun.
Please recommend how to make the nonlinear inequality constraint perform correctly.
Thanks

Best Answer

I am not sure that I understand why you believe that there is a problem with the result. Often, the minimum for a constrained function occurs on the constraint boundaries.
If you think that fmincon is not finding a global minimum, well, look at the suggestions in the documentation for improving results, and look at Searchng for a Smaller Minimum.
Good luck,
Alan Weiss
MATLAB mathematical toolbox documentation
Related Question