Code covered by the BSD License  

Highlights from
OpenOpt

image thumbnail

OpenOpt

by

 

25 Nov 2006 (Updated )

nonSmoothSolve (similar to fsolve), non-smooth & noisy local + some global solvers; works in Octave

nonlinconstr(prob)
function r = nonlinconstr(prob)
global globstat

prob.cnl = [];
prob.cnl.lb = prob.lb;
prob.cnl.ub = prob.ub;

if ~isempty(prob.b) || ~isempty(prob.beq) || ...
        (~isempty(prob.lb) && any(isfinite(prob.lb))) || ...
        (~isempty(prob.ub) && any(isfinite(prob.ub)))
    prob = lin2nonlin(prob);
end

%setting Lagrange multipliers - AFTER lin2nonlin()!

lambda = zeros(prob.nh,1);%beware! - zeros can yield unbounded solution
mu = zeros(prob.nc,1);

maxUpdateMultipliers = 100;
ExtractRoutineParamsFromProb;
if ~isfield(prob, 'lambda'); prob.lambda = lambda; end
if ~isfield(prob, 'mu'); prob.mu = mu; end

FUNCS = {'f' 'd2f' 'c' 'h' 'dc' 'dh' 'd2c' 'd2h'};
% NO df!! - see df handling below!!

prob.cnlf = prob.f;
prob.cnldf = prob.df;
if ~isempty(prob.df) && isempty(prob.primal.df)
    prob.cnl.df = @(x, prob, varargin) prob.cnldf(x, prob, [], prob.cnlf);
    prob.df = @(x, prob, varargin) cnl_df(x, prob, varargin, prob.cnlf);
else
    FUNCS{end+1} = 'df';
end
% for fn = {'df' 'f' 'd2f' 'c' 'h' 'dc' 'dh' 'd2c' 'd2h'} % NO df!! - see df handling above!!
for fn = FUNCS
    if ~isfield(prob, fn{1}) || isempty(getfield(prob, fn{1})); continue; end %#ok<GFLD>
    F = getfield(prob, fn{1}); %#ok<GFLD>
    prob.cnl = setfield(prob.cnl, fn{1}, F); %#ok<SFLD>
    eval(['prob.' fn{1} '= @cnl_' fn{1} ';']); %#ok<SFLD>

    %debug
    if ~ismember(fn, {'f' 'c' 'h' 'df' 'dc' 'dh'}) && ~isempty(getfield(prob.primal, fn{1})) %#ok<GFLD>
        prob.err(['nonlinconstr wrapper: ' fn{1} ' isn''t implemented yet'])
    end
    %debug end
end

%todo: remove (make auto before 1st call)
prob = updateLM(prob.x0, prob);

prob.cnl.iterfcn = prob.iterfcn;
prob.iterfcn = @nonLinConstrIterFcn;
% prob.cnl.iterfcn = prob.iterfcn;
% prob.iterfcn = @LMiter;

prob.nc = 0;
prob.nh = 0;
% todo: add prob.cnl.nc, prob.cnl.nh

nUpdateMultipliers = 0;
prob.f0 = prob.f(prob.x0, prob);
% prob.advanced.ignoreIStop(end+1)=5;
nn=10;
if isequal(lower(prob.solverName), 'hpso')
    prob.MaxIter = prob.MaxIter/nn;
end
for itn=1:maxUpdateMultipliers
    r = prob.solver(prob);
%     if prob.debug; disp(['istop: ' num2str(r.istop)]); end
    %     disp(r.xf)
    if r.istop == 13; return; end
    if ~stopcase(r) % MaxTime, MaxCPUTime, MaxIter, MaxFunEvals or etc is exeeded
        if ~isequal(lower(prob.solverName), 'hpso') || itn==nn
            break
        end
    end

    if ~isempty(prob.c)%todo: if r includes r.c info
        c = prob.cnl.c(r.xf, prob);
    else c = [];
    end
    if ~isempty(prob.h)
        h = prob.cnl.h(r.xf, prob);
    else h = [];
    end
    bViolationH = abs(h)>prob.TolCon;%matrix
    bViolationC = c>prob.TolCon;%matrix
    indViolationH = find(bViolationH);
    indViolationC = find(bViolationC);%#ok<EFIND> %matrix
    if isempty(indViolationH) && isempty(indViolationC) && r.istop ~= -8.8888
        r.advanced.nUpdateMultipliers = nUpdateMultipliers;
        % todo: hint for setting start lambda, mu
        break
    end
    nUpdateMultipliers = nUpdateMultipliers + 1;
    if nUpdateMultipliers>maxUpdateMultipliers
        r.istop = -13;
        return
    end
    %     if prob.advanced.debug
    if prob.iterPrint>0 && ~isequal(lower(prob.solverName), 'hpso')
        prob.info(['nIter: ' int2str(globstat.nIter) '  nUpdateLagrangeMultipliers = ' int2str(nUpdateMultipliers)])
    end

    if prob.debug
        disp(['istop: ' num2str(r.istop)])
        if r.istop == -5; disp('No line-search minimum is found'); end
        s = '';
        if ~isempty(c); s = horzcat(s, '  ', int2str(length(indViolationC)), ' ineq'); end
        if ~isempty(h); s = horzcat(s, '  ', int2str(length(indViolationH)), ' eq'); end
        prob.info(['constraints violation: ' s])
    end

    %     [prob.lambda prob.mu] = updateLM(prob);
    %     if ~isempty(bViolationH)
    %         j = find(bViolationH & ~prob.lambda);
    %         prob.lambda(j) = prob.lambda(j) + 1;
    %     end
    %     if ~isempty(bViolationC)
    %         j = find(bViolationC & ~prob.mu);
    %         prob.mu(j) = prob.mu(j) + 1;
    %     end
%     if r.istop ~= -5
        prob.x0 = r.xf;
%     end
    prob = updateLM(prob.x0, prob);
    %     globstat.xPrev=[];%checkme!
    if r.istop ~= -5%no line-search found
        prob.f0 = prob.f(prob.x0, prob);
    end
    if ismember(prob.solverName, {'ralg' 'hPSO'}) && isfield(r, 'advanced') && isfield(r.advanced, 'ralg') && isfield(r.advanced.ralg, 'hs')
        %         prob.ralg.b = r.advanced.ralg.b;
        prob.ralg.hs = r.advanced.ralg.hs;
    end
    %     prob.f0 = r.ff;
    globstat.isFinished = 0;

    if isfield(globstat, 'lastX')
        globstat.lastX = inf*ones(size(globstat.lastX));
    end

    globstat.ff = inf;
    globstat.bestOptimFVal = inf;
    if isfield(globstat, 'lastOptimFVal')
        globstat.lastOptimFVal = inf;
    end
end
if itn==maxUpdateMultipliers && (~(isempty(indViolationH) && isempty(indViolationC)) || r.istop == -8.8888)
    r.istop = -13;
end
if ~isempty(prob.lambda); r.lambda = prob.lambda; end
if ~isempty(prob.mu); r.mu = prob.mu; end
end

% function prob = updateLM(x, prob)
% %linear constraints must be alreay converted to non-lin
% % f = prob.f(x, prob);
% % ndf = norm(prob.df(x, prob),1);
% DF = prob.cnldf(x, prob);
% % ndf = norm(DF,1);
%
% if isfield(prob.cnl, 'dc')
%     c = prob.cnl.c(x, prob);
% %     Nc = norm(c, 1)+realmin;
%     dc = prob.cnl.dc(x, prob);
% else
%     dc = [];
% end
% if isfield(prob.cnl, 'dh')
%     h = prob.cnl.h(x, prob);
% %     Nh = norm(h, 1)+realmin;
%     dh = prob.cnl.dh(x, prob);
% else
%     dh = [];
% end
%
% p = [c; abs(h)];
% p(p<prob.TolCon) = 0;
% if ~any(p); return; end
% dh2 = dh .* repmat(sign(h), 1, size(dh,1)).';
% A = [dc dh2];
% b = DF;
% % A = sum(A);
% % b = sum(b);
% x0 = [prob.mu; prob.lambda];
% [x lam exitflag output] = linprog(p,A,b,[],[],zeros(length(x0),1),[],[]);
% if exitflag<=0
%     disp(exitflag)
%     disp(output)
%     keyboard
% end
% prob.mu = x(1:length(prob.mu));
% prob.lambda = x(1+length(prob.mu):end);
% end

Contact us