fmincon constraint function not continous

Firstly, based on the Matlab document, 'the constraint functions should be continous and have continuous first deviatives'. My understanding is that the function should be continous over (-/infinite to +/infinite), not over the low bound and up bound that I input in the code. For example, my lb=1 and ub=100 and over this interval, the function y=1/x is continous, but I still cant directly use it.

Seconldy, if my understanding above is right, I only have one choice to transfer the constraint function to continous, take the previous example again, I need to input x*y=1 instead of y=1/x even though my low bound is 1. But the problem is that what if my constraint function is complicated and I cant transfer it to a continous form. Currently I am trying to minimize a problem stated below, and my understanding is that even though my constraint function is not continous, I can use many different initial point to address this issue (i.e. use different initial point to get the optimised results and then that the minimum one), I dont know whether my understanding is right. My matlab code is as follows (I use 60 different initial points).

load('initial_points.mat','matrix')
initial_point=matrix; % this is a 60 by 2 matrix

D_intial_total=zeros(1,5); %this is the matrix to keep record of all the results (not to be used in this code).

%%to get different results for different 'a' value. (we use different intial points to get local minimise results and then further take the one with the minimum costs.)

k=1;
for a=1:20
    fun = @(x)a*x(1)+x(2);
    lb = [0.01,0.01];
    ub = [1/9-0.00001,4];
    A = [];
    b = [];
    Aeq = [];
    beq = [];
    ii=1;
    for initialindex=1:size(initial_point,1)
        x0 = initial_point(initialindex,:);
        options = optimoptions('fmincon','Algorithm','sqp');
        %'interior-point'  'trust-region-reflective''sqp''sqp-legacy''active-set'
        % This is the constrate of \zeta and u_t where x(1) is u_t and x(2) is
        % u_v
        c = @(x)real((-2 + 2*x(1) - 3*x(2) + 3*x(1)*x(2) - 3*x(2)^2)/(3*x(2)^3) - ((1 + sqrt(3)*i)*(-3*x(2)^3*(1 - x(1) - 3*x(1)*x(2)) - (-2 + 2*x(1) - 3*x(2) +3*x(1)*x(2) - 3*x(2)^2)^2))/(3*2^(2/3)*x(2)^3*(16 - 48*x(1) + 48*x(1)^2 - 16*x(1)^3 + 72*x(2) - 216*x(1)*x(2) + 216*x(1)^2*x(2) - 72*x(1)^3*x(2) +180*x(2)^2 - 468*x(1)*x(2)^2 + 396*x(1)^2*x(2)^2 - 108*x(1)^3*x(2)^2 + 288*x(2)^3 - 630*x(1)*x(2)^3 + 396*x(1)^2*x(2)^3 - 54*x(1)^3*x(2)^3 + 297*x(2)^4 - 540*x(1)*x(2)^4 + 243*x(1)^2*x(2)^4 + 189*x(2)^5 - 270*x(1)*x(2)^5 + 81*x(1)^2*x(2)^5 + 54*x(2)^6 - 54*x(1)*x(2)^6 + sqrt((16 - 48*x(1) + 48*x(1)^2 -16*x(1)^3 + 72*x(2) - 216*x(1)*x(2) + 216*x(1)^2*x(2) - 72*x(1)^3*x(2) + 180*x(2)^2 - 468*x(1)*x(2)^2 + 396*x(1)^2*x(2)^2 - 108*x(1)^3*x(2)^2 + 288*x(2)^3-630*x(1)*x(2)^3 + 396*x(1)^2*x(2)^3 - 54*x(1)^3*x(2)^3 + 297*x(2)^4 - 540*x(1)*x(2)^4 + 243*x(1)^2*x(2)^4 + 189*x(2)^5 - 270*x(1)*x(2)^5 + 81*x(1)^2*x(2)^5+54*x(2)^6 - 54*x(1)*x(2)^6)^2 + 4*(-3*x(2)^3*(1 - x(1) - 3*x(1)*x(2)) - (-2 + 2*x(1) - 3*x(2) + 3*x(1)*x(2) - 3*x(2)^2)^2)^3))^(1/3)) +((1 - sqrt(3)*i)*(16 - 48*x(1) + 48*x(1)^2 - 16*x(1)^3 + 72*x(2) - 216*x(1)*x(2) + 216*x(1)^2*x(2) - 72*x(1)^3*x(2) + 180*x(2)^2 - 468*x(1)*x(2)^2 + 396*x(1)^2*x(2)^2 -108*x(1)^3*x(2)^2 + 288*x(2)^3 - 630*x(1)*x(2)^3 + 396*x(1)^2*x(2)^3 - 54*x(1)^3*x(2)^3 + 297*x(2)^4 - 540*x(1)*x(2)^4 + 243*x(1)^2*x(2)^4 + 189*x(2)^5 - 270*x(1)*x(2)^5 +81*x(1)^2*x(2)^5 + 54*x(2)^6 - 54*x(1)*x(2)^6 + sqrt((16 - 48*x(1) + 48*x(1)^2 - 16*x(1)^3 + 72*x(2) - 216*x(1)*x(2) + 216*x(1)^2*x(2) - 72*x(1)^3*x(2) + 180*x(2)^2 -468*x(1)*x(2)^2 + 396*x(1)^2*x(2)^2 - 108*x(1)^3*x(2)^2 + 288*x(2)^3 - 630*x(1)*x(2)^3 + 396*x(1)^2*x(2)^3 - 54*x(1)^3*x(2)^3 + 297*x(2)^4 - 540*x(1)*x(2)^4 +243*x(1)^2*x(2)^4 + 189*x(2)^5 - 270*x(1)*x(2)^5 + 81*x(1)^2*x(2)^5 + 54*x(2)^6 - 54*x(1)*x(2)^6)^2 + 4*(-3*x(2)^3*(1 - x(1) - 3*x(1)*x(2)) - (-2 + 2*x(1) -3*x(2) + 3*x(1)*x(2) - 3*x(2)^2)^2)^3))^(1/3))/(6*2^(1/3)*x(2)^3))-0.1+1e-6; 
        ceq = [];
        nonlcon = @(x)deal(c(x),[]);
        x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon,options);
        % use D_intial_point_result to storage the optimisation results for one
        % 'a' value with different intial points, then further find the optimised
        % results from 'D_intial_point_result'.
        D_intial_point_result(ii,5)=a;
        D_intial_point_result(ii,3)=a*x(1)+x(2); % this is the cost
        D_intial_point_result(ii,1:2)=x; 
        D_intial_point_result(ii,4)=real((-2 + 2*x(1) - 3*x(2) + 3*x(1)*x(2) - 3*x(2)^2)/(3*x(2)^3) - ((1 + sqrt(3)*i)*(-3*x(2)^3*(1 - x(1) - 3*x(1)*x(2)) - (-2 + 2*x(1) - 3*x(2) +3*x(1)*x(2) - 3*x(2)^2)^2))/(3*2^(2/3)*x(2)^3*(16 - 48*x(1) + 48*x(1)^2 - 16*x(1)^3 + 72*x(2) - 216*x(1)*x(2) + 216*x(1)^2*x(2) - 72*x(1)^3*x(2) +180*x(2)^2 - 468*x(1)*x(2)^2 + 396*x(1)^2*x(2)^2 - 108*x(1)^3*x(2)^2 + 288*x(2)^3 - 630*x(1)*x(2)^3 + 396*x(1)^2*x(2)^3 - 54*x(1)^3*x(2)^3 + 297*x(2)^4 - 540*x(1)*x(2)^4 + 243*x(1)^2*x(2)^4 + 189*x(2)^5 - 270*x(1)*x(2)^5 + 81*x(1)^2*x(2)^5 + 54*x(2)^6 - 54*x(1)*x(2)^6 + sqrt((16 - 48*x(1) + 48*x(1)^2 -16*x(1)^3 + 72*x(2) - 216*x(1)*x(2) + 216*x(1)^2*x(2) - 72*x(1)^3*x(2) + 180*x(2)^2 - 468*x(1)*x(2)^2 + 396*x(1)^2*x(2)^2 - 108*x(1)^3*x(2)^2 + 288*x(2)^3-630*x(1)*x(2)^3 + 396*x(1)^2*x(2)^3 - 54*x(1)^3*x(2)^3 + 297*x(2)^4 - 540*x(1)*x(2)^4 + 243*x(1)^2*x(2)^4 + 189*x(2)^5 - 270*x(1)*x(2)^5 + 81*x(1)^2*x(2)^5+54*x(2)^6 - 54*x(1)*x(2)^6)^2 + 4*(-3*x(2)^3*(1 - x(1) - 3*x(1)*x(2)) - (-2 + 2*x(1) - 3*x(2) + 3*x(1)*x(2) - 3*x(2)^2)^2)^3))^(1/3)) +((1 - sqrt(3)*i)*(16 - 48*x(1) + 48*x(1)^2 - 16*x(1)^3 + 72*x(2) - 216*x(1)*x(2) + 216*x(1)^2*x(2) - 72*x(1)^3*x(2) + 180*x(2)^2 - 468*x(1)*x(2)^2 + 396*x(1)^2*x(2)^2 -108*x(1)^3*x(2)^2 + 288*x(2)^3 - 630*x(1)*x(2)^3 + 396*x(1)^2*x(2)^3 - 54*x(1)^3*x(2)^3 + 297*x(2)^4 - 540*x(1)*x(2)^4 + 243*x(1)^2*x(2)^4 + 189*x(2)^5 - 270*x(1)*x(2)^5 +81*x(1)^2*x(2)^5 + 54*x(2)^6 - 54*x(1)*x(2)^6 + sqrt((16 - 48*x(1) + 48*x(1)^2 - 16*x(1)^3 + 72*x(2) - 216*x(1)*x(2) + 216*x(1)^2*x(2) - 72*x(1)^3*x(2) + 180*x(2)^2 -468*x(1)*x(2)^2 + 396*x(1)^2*x(2)^2 - 108*x(1)^3*x(2)^2 + 288*x(2)^3 - 630*x(1)*x(2)^3 + 396*x(1)^2*x(2)^3 - 54*x(1)^3*x(2)^3 + 297*x(2)^4 - 540*x(1)*x(2)^4 +243*x(1)^2*x(2)^4 + 189*x(2)^5 - 270*x(1)*x(2)^5 + 81*x(1)^2*x(2)^5 + 54*x(2)^6 - 54*x(1)*x(2)^6)^2 + 4*(-3*x(2)^3*(1 - x(1) - 3*x(1)*x(2)) - (-2 + 2*x(1) -3*x(2) + 3*x(1)*x(2) - 3*x(2)^2)^2)^3))^(1/3))/(6*2^(1/3)*x(2)^3))-0.1;   
        % to check whether the optimised result satisfy constrait
        ii=ii+1;
    end
    D_intial_total=[D_intial_total;D_intial_point_result];  % to keep a record of all optimised reults (not to be used in this code)
    D_intial_point_result=D_intial_point_result(find(D_intial_point_result(:,4)<0),:); % remove the optimised result that is not satisfy the constrait 
    [min_value, min_index] = min(D_intial_point_result(:,3)); % find global minimiser
    % opresultutlinear is the main result
    opresultutlinear(k,1)=a;  
    opresultutlinear(k,2:3)=D_intial_point_result(min_index,1:2);% u_t and \zeta 
    opresultutlinear(k,4)=D_intial_point_result(min_index,4); % the constrait value
    k=k+1;
end

Answer 1

[Walter Roberson]

My understanding is that the function should be continous over (-/infinite to +/infinite), not over the low bound and up bound that I input in the code.

No, it has to be continuous within your bounds. Gradients will be estimated numerically, which can validly be done if the function is continuous with continuous first derivatives everywhere within the range being tested.

[Tianshu Gao]

Thank you for your answer.

What you said really make sense to me since the algorithm should only work winthin the bounds. But the discription in Matlab document worries me (it said that 'fmincon is a gradient-based method that is designed to work on problems where the objective and constraint functions are both continuous and have continuous first derivatives.'). Do you think it is just because they made mistakes when they write this discription (they forget to mention that continuous only need to within the range) ?

[Torsten]

They didn't mention it because for people who are concerned with optimization, it's self-evident.

By the way: The function y = 1/x is continuous everywhere where it is defined. x = 0 is not a discontinuity, but a point where the function is undefined.

[Matt J]

x = 0 is not a discontinuity, but a point where the function is undefined.

That would also make it discontinuous at x=0.

[Matt J]

No, it has to be continuous within your bounds.

Unless, you are using some non-default setting which allows the iterations to stray outside the bounded region, like the interior-point method with HonorBounds=false.

[Torsten]

Ok, it's hairsplitting, but in a mathematical sense, the terms "continuous" and "discontinuous" only make sense for points where the function is defined. 0 in this case is a singular point.

[Tianshu Gao]

Thank you for all above comments. Now I understand the my question 1. But how about my question 2, if I have a very complited expression like the one in the code, what if I cant prove it is continuous or not.

[Torsten]

If you don't have jumps in the coefficients of the expression, if you have no if-statements that determine the expression for different ranges, if there are no integer solution variables in the expression I'd say: just assume it is continuous and see if the solver returns something useful.