bandb.m

数学优化工具箱

function [optimal,itercount]=bandb(opts)
    % main function for the l1l1 branch and bound algorithm
    % this takes a structure with all the problem specs as input
    % and returns a structure containing the optimal output

    optimal=newcube(opts.e);
    candidates = cell(1);

    n = opts.n;
    itercount = 0;

    % bulid the initial sedumi structures and rectangle
    cube0 = bandb_initial_cube(opts);
    cube0 = bound(cube0,opts);

    if cube0.feasible == 0;
        fprintf('L1L1:The starting region is infeasible\n');
        fprintf('Exiting\n');
    else
        fprintf('L1L1:The starting region is feasible\n');
        candidates{1}=cube0;
        optimal = cube0;

        while 1;
            newcandidates = cell(0);
            minlb = 2*opts.e;
            mini = -1;
            gap = 1.0;


            lbgap = 1e10;
       
            for i = 1:length(candidates);

                lbgap = min(lbgap,candidates{i}.lb);
               % fprintf('%f %f %f %f %f\n',candidates{i}.lb,optimal.v,opts.epsilon*optimal.v,candidates{i}.lb/optimal.v,opts.epsilon);
                if candidates{i}.lb < opts.epsilon*optimal.v;
                    newcandidates{end+1}=candidates{i};
                    if candidates{i}.lb < minlb;
                        mini = length(newcandidates);
                        minlb = candidates{i}.lb;

                    end
                end
            end
            %    fprintf('lb gap %f\n',optimal.v-lbgap);
            %printf('Length : %d\n',numel(newcandidates));
            if (length(newcandidates) == 0) | (itercount > opts.maxiter);

                if (optimal.v > opts.e);
                    %    fprintf('Optimal Solution found\n');
                    %else
                    fprintf('No feasible point found\nExiting..\n');
                end
                break;
            end

            ol = length(candidates);

            bestCube =  newcandidates{mini};
            candidates = newcandidates;
            candidates(mini)=[];
            nl= length(candidates);

            [cubel,cuber] = branch(bestCube,opts);

            if cubel.feasible ==1;
                cubel = bound(cubel,opts);
                if cubel.feasible == 1;
                    candidates{end+1}=cubel;
                    if (cubel.v < optimal.v);
                        optimal=cubel;
                    end
                end
            end;

            if cuber.feasible ==1;
                cuber = bound(cuber,opts);
                if cuber.feasible == 1;
                    candidates{end+1}=cuber;
                    if (cuber.v < optimal.v);
                        optimal=cuber;
                    end
                end
            end

            lbgap = min([lbgap,cubel.lb,cuber.lb]);

            gap = lbgap/optimal.v;

            fprintf('[% 3d:% 3d][%8f] %8f: (%8f %8f %8f) | (%8f %8f %8f)\n',itercount,length(candidates),gap,optimal.v, cubel.lb,cubel.v,cubel.obj,cuber.lb,cuber.v,cuber.obj);
            %fprintf('[% 3d:% 3d] %8f(%8f)\n',itercount,length(candidates),optimal.v,gap);
            itercount=itercount+1;
        end

        gap = lbgap/optimal.v;
        fprintf('Final:    %8f(%8f)\n',optimal.v,gap);

    end
return

function cube = bandb_initial_cube(opts);
		cube=newcube(opts.e);
		cube.H=opts.ybounds;
return

function [cubel,cuber]=branch(cube,opts);
	% perform the branching operation for the bandb algorithm.

	alpha = opts.alpha;
	cubel=newcube(opts.e);
	cuber=newcube(opts.e);
	H = cube.H;
	cubel.H=H;
	cuber.H=H;

	if opts.branch_type ==1;
	  width = (H(:,2)-H(:,1)).*opts.branch_array;
	else
	  width = (H(:,2)-H(:,1));
	end

	[maxw,maxi]=max(width);
	newborder= H(maxi,1)+alpha*maxw;

	cubel.H(maxi,2)=newborder;
	cubel.bdim = maxi;

	cubel = bandb_refine_bounds(cubel,opts);
	cubel.refined = 1;

	cuber.H(maxi,1)=newborder;
	cuber.bdim = maxi;
	cuber = bandb_refine_bounds(cuber,opts);
	cuber.refined = 1;	
return

function value = newcube(e);

if nargin < 1;
    e = 100000000;
end

value ={};
value.H=[];
value.point=[];
value.feasible=0;
value.refined = 0;

value.lbr=[];value.lb = e;
value.objr=[];value.obj=e+1;
value.vr=[];value.v=e+2;
value.den=[];
value.bdim = 0;
return