global search with fmincon returns [] when a particular constraint is added

Hello, I'm running a global search algorithm to determine the constraint optimum of a loglikelihood function. I'm repeating the same optimization for all k < R, in order to populate the matrix xfinals

    % Create the control matrix X with all the required control variables:
    jj = ones((C-2),1);                         % vector of constants
    Xscale = [ jj,   matrix(:,4),matrix(:,5),matrix(:,6),matrix(:,7),matrix(:,8),matrix(:,9),matrix(:,11),matrix(:,12),matrix(:,13),matrix(:,14),matrix(:,15),matrix(:,16),matrix(:,18), matrix(:,19),matrix(:,20), matrix(:,23), matrix(:,24), matrix(:,25), matrix(:,26) ];    % Fixed effect for scale parameter lambda
Xshape = [ jj, matrix(:,4),matrix(:,15),matrix(:,16),matrix(:,23), matrix(:,24), matrix(:,25), matrix(:,27) ];                  % Fixed effect for shape parameter k

    % Fixed effect for scale parameter lambda
    %Xscale = [ jj, matrix(:,4), matrix(:,5),matrix(:,7), matrix(:,9), matrix(:,11), matrix(:,12), matrix(:,13), matrix(:,14), matrix(:,15),matrix(:,16),matrix(:,18), matrix(:,19), matrix(:,20), matrix(:,21),matrix(:,24), matrix(:,23), matrix(:,26),matrix(:,27) ];
    % Fixed effect for shape parameter k
    %Xshape = [ jj, matrix(:,14), matrix(:,19), matrix(:,20), matrix(:,21),matrix(:,24), matrix(:,23), matrix(:,26),matrix(:,27) ];

    [RXscale, CXscale] = size(Xscale);
    [RXshape, CXshape] = size(Xshape);

    % Total number of unknown per alpha, beta combination:
    unknowns = CXscale + CXshape + 1;

    % The initial points (theta0) is a matrix with 3 rows (one per parameter)
    % and a number of columns equal to the number the number of variables in X.
    theta0 = rand(unknowns,1) * .00005;%*.000001;         % initial points (column vector) * 0.001

    Ae = [];
    be = [];
    Aeq = [];
    beq = [];
    lb = - ones(1, unknowns) * 100;
    ub = ones(1, unknowns) * 100;
    for k=1:R
    if count_min_v(k,1) > 80
        xfinals(k,:) = zeros(1,unknowns);
        loglikelihood(k,1) = 0;
        exitreason(k,1) = 6;
    else
        % make sure to transpose the vector of bids to be sent to the MLE
        % function

        % Set MaxIter and MaxfunEvals to 1000
        values = (v_pos(k,1:(C-2)))';
        f = @(theta) fun_MLE_controls( values, theta, Xscale, CXscale, Xshape, CXshape );
        opts = optimoptions(@fmincon,'Display','final-detailed','MaxIter', 20000,...
            'MaxfunEvals',20000,'TolX',1e-18,'TolFun',1e-18,'TolCon',1e-18,'Algorithm','interior-point','Hessian','bfgs');
        mycon = @(theta) constraint( values, theta, Xscale, CXscale, Xshape, CXshape ) ;
        problem = createOptimProblem('fmincon','objective',f,'x0',theta0,'lb',lb,'ub',ub,'nonlcon',mycon,'options',opts) ;

        gs = GlobalSearch('Display','iter','TolFun',1e-18,'TolX',1e-18,'StartPointsToRun','bounds-ineqs');
        [xfinal, fval, exitflag, output] = run(gs,problem);
        %[xfinal, fval, exitflag] = fmincon(f,theta0,Ae,be,Aeq,beq,lb,ub,mycon,options);        % no nonlinear constraint []
        xfinals(k,:) = xfinal;            
        exitreason(k,1) = exitflag;
        loglikelihood(k,1) = fval; %fun_MLE_controls( values, xfinals(k,:), Xscale, CXscale, Xshape, CXshape );

        end
    end

Since I want the expression in c(2) and c(3) to be larger than zero for all values, I tell the code to make the smallest value of the vector operation larger than zero. If the smallest value is larger than zero, then all the values in the vector are lareger than zero.

    function [ c,ceq ] = constraint( v, parameters, Xscale, CXscale, Xshape, CXshape )

    % LAMBDA and K as in fun_MLE_controls must be larger than zero

    % Define the Gamma-Weibull parameters
    lambda = parameters(1:CXscale,1);               % scale
    k = parameters((CXscale+1):(CXshape+CXscale),1);        % shape        
    zi = parameters((CXscale+CXshape+1),1);             % unobserved heterogeneity

    LAMBDA = exp(sum(Xscale * lambda));
    K = exp(sum(Xshape * k));

    lk = v.^K;              % lk = v_i^k
    sum_k = sum(lk);

    THETA = LAMBDA^(-K);
    c(1) = - ((K + zi)/K);
    pdf = K * THETA^(1+(zi/K)) * (v.^(zi+K-1)) .* exp(-THETA * v.^K) / gamma(1+(zi/K));
    c(2) = -min(pdf);
    c(3) = -min((v.^(-1)).*((zi + K - 1) - K .* THETA .* (v.^K)) + ...
    (1./(gammainc(THETA .* v.^K,(1 + (zi/K)),'upper'))));
    ceq = [];

I do not report the MLE function, but there is no problem with that because the code works fine without constraint c(3). When I include c(3) xfinal = [] and the program stops with the following message:

    Subscripted assignment dimension mismatch.

    Error in fioretti_auction_Charity_controls_GlobalSearch (line 270)
        xfinals(k,:) = xfinal;

Could anyone help me please?

Answer 1

[Matt J]

You need to decide what goes in xfinals when no feasible solution is found. Something like this would make sense:

    if ~isempty(xfinal)
     xfinals(k,:) = xfinal;     

    else
      xfinals(k,:)=nan;
    end

[Mic]

Thank you Matt for your input, but I do not understand why it returns []. It should at least return the starting points right? Also, the problem is in the "gammainc" function, because if I remove that function the code goes through without problem.

[Matt J]

Thank you Matt for your input, but I do not understand why it returns []. It should at least return the starting points right?

What exitflag value is returned in the cases where xfinal=[]? That will help us interpret what is being returned for xfinal.

From the documentation on GlobalSearch, it doesn't look to me that a solution will be accepted if fmincon does not return a positive exitflag, indicating convergence. Thus, it stands to reason that if no point satisfying the constraints exists, the set of solutions returned by run() could be empty.

Also, the problem is in the "gammainc" function, because if I remove that function the code goes through without problem.

That suggests that you are over-constraining your problem, i.e., the addition of the gammainc constraint makes the whole problem infeasible. But, we need the exitflag to know what's really going on.

Another thing I now notice is that your tolerances TolCon, TolFun, etc... have been set very tight (1e-18). That could make it very difficult to achieve the convergence criteria, which has the same effect as if the problem is infeasible.