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

updateLM(x, prob)
function prob = updateLM(x, prob)
% linear constraints must be alreay converted to non-lin
% f = prob.f(x, prob);
normcase = 1;
power = 2;
reduce = 0.75;
ExtractRoutineParamsFromProb;

if isfield(prob.cnl, 'dc')
    c = prob.cnl.c(x, prob);
    dc = prob.cnl.dc(x, prob);
    DC = dc(:,c>prob.TolCon);
    MaxC = max(c);
else
    [c dc DC MaxC] = deal([]);
end
if isfield(prob.cnl, 'dh')
    h = abs(prob.cnl.h(x, prob));
    dh = prob.cnl.dh(x, prob);
    DH = dh(:,h>prob.TolCon);
    MaxH = max(h);
else
    [h dh DH MaxH] = deal([]);
end

MaxConstr = max([MaxC MaxH]);
ndf = norm(prob.cnldf(x, prob),normcase);
%keyboard
if ~ndf % for example, sometimes from nonSmoothSolve()
    NDHm = norm(DH(:),normcase);
    NDCm = norm(DC(:),normcase);%avoid zerodevide
NDAll = [NDCm NDHm];

%fixing for Octave (to not produce warn "division by zero")
if isempty(NDAll); ndf = 0; else ndf = mean(NDAll); end

%if isempty(NDHm); NDHm = 0; end
%if isempty(NDCm); NDCm = 0; end
%    ndf = mean(NDCm + NDHm)/(size(DC,2)+size(DH,2)+realmin); 
end 
%keyboard

switch prob.advanced.cPenalty
    case 'linear'
        for i = 1:size(dc,2)
            if c(i) <= prob.TolCon; continue; end
            ndc = norm(dc(:,i),normcase);
            if ~ndc; continue; end
            if length(h)+length(c)>1
                K = (2.05-(1-c(i)/MaxConstr)^power);
            else
                K=1.05;
            end
            prob.mu(i) = max(reduce * prob.mu(i) + (1-reduce)*K * ndf./ndc, K * ndf./ndc);
        end
    case 'norm2'
%         prob.mu = max(ndf./NDCmean, norm(prob.mu,1)*1.1);
    otherwise
        prob.err('undefined or unimplemented cPenalty type')
end
switch prob.advanced.hPenalty
    case 'linear'
        for i = 1:size(dh,2)
            if abs(h(i)) <= prob.TolCon; continue; end
            ndh = norm(dh(:,i),normcase);
            if ~ndh; continue; end
                if length(h)+length(c)>1
                    K = (2.05-(1-h(i)/MaxConstr)^power);
                else
                    K=1.05;
                end
                prob.lambda(i) = max(reduce*prob.lambda(i) + (1-reduce) * ndf./ndh * K, ndf./ndh * K);
        end
    case 'norm2'
%         prob.lambda = max(ndf./NDHmean, norm(prob.lambda,1)*1.5);
    otherwise
        prob.err('undefined or unimplemented hPenalty type')
end
end

Contact us