addevalvariablecuts.m

matlab波形优化算法经常要用到的matlab toolbox工具箱:yalmip

function pcut = addEvalVariableCuts(p)

temp = p.F_struc(1:p.K.f,:);
p.F_struc(1:p.K.f,:) = [];
pcut = p;
for i = 1:length(p.evalMap)
    y = p.evalVariables(i);
    x = p.evalMap{i}.variableIndex;
    xL = p.lb(x);
    xU = p.ub(x);

    % Generate a convex hull polytope
    if xL<xU
        if ~isempty(p.evalMap{i}.properties.convexhull)
            % A convex hull generator function is available!
            K = [];
            try
                [Ax,Ay,b,K]=feval(p.evalMap{i}.properties.convexhull,xL,xU, p.evalMap{i}.arg{2:end-1});
            catch
                [Ax,Ay,b]=feval(p.evalMap{i}.properties.convexhull,xL,xU, p.evalMap{i}.arg{2:end-1});
            end
        else
            if length(xL)>1
                disp(['The ' p.evalMap{i}.fcn ' operator does not have a convex hull operator'])
                disp('This is required for multi-input single output operators');
                disp('Sampling approximation does not work in this case.');
                error('Missing convex hull operator');
            end
            % sample function
            z = linspace(xL,xU,100);
            
            if isequal(p.evalMap{i}.fcn,'power_internal2')
                % Special code for automatically converting sigmonial
                % terms to be solvable with bmibnb
                fz = feval(p.evalMap{i}.fcn,z,p.evalMap{i}.arg{2});

            else
                arg = p.evalMap{i}.arg;
                arg{1} = z;
                fz = feval(p.evalMap{i}.fcn,arg{1:end-1});
          % end
            [minval,minpos] = min(fz);
            [maxval,maxpos] = max(fz);
            xtestmin = linspace(z(max([1 minpos-5])),z(min([100 minpos+5])),100);
            xtestmax = linspace(z(max([1 maxpos-5])),z(min([100 maxpos+5])),100);
            arg{1} = xtestmin;
            fz1 = feval(p.evalMap{i}.fcn,arg{1:end-1});
            arg{1} = xtestmax;
            fz2 = feval(p.evalMap{i}.fcn,arg{1:end-1});
            z = [z(:);xtestmin(:);xtestmax(:)];
            fz = [fz(:);fz1(:);fz2(:)];
            [z,sorter] = sort(z);
            fz = fz(sorter);
           [z,ii,jj]=unique(z);
           fz = fz(ii);
            end

            %                    fz = feval(p.evalMap{i}.fcn,z);
            %   end
            % create 4 bounding planes
            % f(z) < k1*(x-XL) + f(xL)
            % f(z) > k2*(x-XL) + f(xL)
            % f(z) < k3*(x-XU) + f(xU)
            % f(z) > k4*(x-XU) + f(xU)
            k1 = max((fz(2:end)-fz(1))./(z(2:end)-xL))+1e-12;
            k2 = min((fz(2:end)-fz(1))./(z(2:end)-xL))-1e-12;
            k3 = min((fz(1:end-1)-fz(end))./(z(1:end-1)-xU))+1e-12;
            k4 = max((fz(1:end-1)-fz(end))./(z(1:end-1)-xU))-1e-12;
            Ax = [-k1;k2;-k3;k4];
            Ay = [1;-1;1;-1];
            b =  [k1*(-z(1)) + fz(1);-(k2*(-z(1)) + fz(1));k3*(-z(end)) + fz(end);-(k4*(-z(end)) + fz(end))];
            
        end
        if ~isempty(b)
            if isempty(K)
                % Compatibility with old code
                K.f = 0;
                K.l = length(b);
            end                
            F_structemp = zeros(size(b,1),length(p.c)+1);
            F_structemp(:,1+y) = -Ay;
            F_structemp(:,1+x) = -Ax;
            F_structemp(:,1) = b;
            pcut.F_struc = [F_structemp; pcut.F_struc];
            pcut.K.l = pcut.K.l + K.l;%length(b);
            pcut.K.f = pcut.K.f + K.f;%length(b);
        end
    end
end
pcut.F_struc = [temp;pcut.F_struc];