mtimes.m

国外专家做的求解LMI鲁棒控制的工具箱,可以相对高效的解决LMI问题

                        % Ok, braze your self for some horrible special case
                        % treatment etc...
                        if isemptynew_mt_hash | bilinearproduct % only search among old monomials
                            if bilinearproduct & candofastlocation
                                before = allbefore_in_old(acounter);
                                if before==0
                                    before = [];
                                else
                                    before = possibleOld(before);
                                end
                            else
                                before = possibleOld(findhash(possibleOldHash,current_hash));                                
                            end
                        else
                            before = findhash(new_mt_hash,current_hash); % first among new monomials
                            if before%~isempty(before) % Try new if we failed finding any match
                                before=before+size(mt_hash,1);
                            else
                                before = possibleOld(findhash(possibleOldHash,current_hash));
                                %                            else
                                % before=before+size(mt_hash,1);
                            end
                        end
                        if before%~isempty(before)
                            Z.lmi_variables(i) = before;
                        else
                            changed_mt=1;
                            isemptynew_mt_hash=0;
                            thesewhereactuallyused(acounter) = 1;
                            if ~bilinearproduct % We don't need to update until later since it not is used inside the loop
                                new_mt_hash = [new_mt_hash;current_hash];
                            end
                            nvar = nvar + 1;
                            Z.lmi_variables(i) = nvar;
                            %  else
                            % Z.lmi_variables(i) = before;
                        end
                        i = i+1;
                    end
                end
                if all(copytofrom)
                    Z.basis = [Z.basis allprodbase];
                else
                    Z.basis = [Z.basis allprodbase(:,find(copytofrom))];
                end
                if all(thesewhereactuallyused)
                    new_mt = [new_mt allmt_xplusy];
                else
                    new_mt = [new_mt allmt_xplusy(:,find(thesewhereactuallyused))];
                end
                if bilinearproduct
                    new_mt_hash = [new_mt_hash;allhash(find(thesewhereactuallyused))];
                end
            end

            if ~isempty(new_mt)
                [i1,j1,k1] = find(mt);
                [ii1,jj1,kk1] = find(new_mt');
                uv = find(used_variables);uv=uv(:);
                mt = sparse([i1(:);ii1(:)+size(mt,1)],[j1(:);uv(jj1(:))],[k1(:);kk1(:)],size(mt,1)+size(new_mt,2),size(mt,2));
            end

            % We pre-allocated a sufficiently long, now pick the ones we
            % actually filled with values
            Z.lmi_variables = Z.lmi_variables(1:i-1);

            % Fucked up order (lmi_variables should be sorted)
            if any(diff(Z.lmi_variables)<0)
                [i,j]=sort(Z.lmi_variables);
                Z.basis = [Z.basis(:,1) Z.basis(:,j+1)];
                Z.lmi_variables = Z.lmi_variables(j);
            end

            [un_Z_vars2] = uniquestripped(Z.lmi_variables);
            if length(un_Z_vars2) < length(Z.lmi_variables)
                [un_Z_vars,hh,jj] = unique(Z.lmi_variables);
                if length(Z.lmi_variables) ~=length(un_Z_vars)
                    Z.basis = Z.basis*sparse([1 1+jj],[1 1+(1:length(jj))],ones(1,1+length(jj)))';
                    Z.lmi_variables = un_Z_vars;
                end
            end

            if changed_mt%~isequal(mt,oldmt)
                newmt = mt(size(oldmt,1)+1:end,:);
                nonlinear = ~(sum(newmt,2)==1 & sum(newmt~=0,2)==1);
                newvariabletype = spalloc(size(newmt,1),1,nnz(nonlinear))';
                nonlinearvariables = find(nonlinear);
                newvariabletype = sparse(nonlinearvariables,ones(length(nonlinearvariables),1),3,size(newmt,1),1)';
                if ~isempty(nonlinear)
                    %mt = internal_sdpvarstate.monomtable;
                    %newvariabletype(nonlinear) = 3;
                    quadratic = sum(newmt,2)==2;
                    newvariabletype(quadratic) = 2;
                    bilinear = max(newmt,[],2)<=1;
                    newvariabletype(bilinear & quadratic) = 1;
                    sigmonial = any(0>newmt,2) | any(newmt-fix(newmt),2);
                    newvariabletype(sigmonial) = 4;
                end
                yalmip('setmonomtable',mt,[oldvariabletype newvariabletype],[mt_hash;new_mt_hash],oldhash);
            end

            if ~(x_isscalar | y_isscalar)
                Z.n = X.n;
                Z.m = Y.m;
            else
                Z.n = max(X.n,Y.n);
                Z.m = max(X.m,Y.m);
            end
        catch
            error(lasterr)
        end
        Z = clean(Z);

    case 2

        n_X = X.n;
        m_X = X.m;
        [n_Y,m_Y] = size(Y);

        x_isscalar =  (n_X*m_X==1);
        y_isscalar =  (n_Y*m_Y==1);

        if ~x_isscalar
            if ((m_X~= n_Y & ~y_isscalar))
                error('Inner matrix dimensions must agree.')
            end
        end

        n = n_X;
        m = m_Y;
        Z = X;

        if x_isscalar
            Z.n = n_Y;
            Z.m = m_Y;
            Z.basis = kron(Z.basis,Y(:));
            Z = clean(Z);
            return
        elseif y_isscalar
            Z.n = n_X;
            Z.m = m_X;
            Z.basis = Z.basis*Y;
            Z = clean(Z);
            return
        end

        Z.n = n;
        Z.m = m;
        Z.basis = kron(Y.',speye(n_X))*X.basis;
        Z = clean(Z);

    case 1

        n_Y = Y.n;
        m_Y = Y.m;
        [n_X,m_X] = size(X);

        x_isscalar =  (n_X*m_X==1);
        y_isscalar =  (n_Y*m_Y==1);

        if ~x_isscalar
            if ((m_X~= n_Y & ~y_isscalar))
                error('Inner matrix dimensions must agree.')
            end
        end

        n = n_X;
        m = m_Y;
        Z = Y;

        % Special cases
        if x_isscalar
            Z.n = n_Y;
            Z.m = m_Y;
            try
                Z.basis = sparse(X)*Y.basis;
            catch
                % This works better when low on memory in some cases
                [i,j,k] = find(Y.basis);
                Z.basis = sparse(i,j,X*k,size(Y.basis,1),size(Y.basis,2));
            end

            Z = clean(Z);
            return
        elseif y_isscalar
            Z.n = n_X;
            Z.m = m_X;
            Z.basis = X(:)*Y.basis;
            Z = clean(Z);
            return
        end

        if m_Y==1
            Z.basis = X*Y.basis;
        else
            try
                speyemy = speye(m_Y);
                kronX = kron(speyemy,X);
                Z.basis = kronX*Y.basis;
            catch
                for i = 1:size(Y.basis,2);
                    dummy = X*reshape(Y.basis(:,i),Y.n,Y.m);
                    Z.basis(:,i) = dummy(:);
                end
            end
        end
        Z.n = n;
        Z.m = m;
        Z = clean(Z);

    otherwise
        error('Logical error in mtimes. Report bug')
end



function Z=clean(X)
temp = any(X.basis,1);
temp = temp(2:end);
index = find(temp);
if ~isempty(index)
    Z = X;
    if length(index)~=length(Z.lmi_variables)
        Z.basis = Z.basis(:,[1 1+index]);
        Z.lmi_variables = Z.lmi_variables(index);
    end
else
    Z = full(reshape(X.basis(:,1),X.n,X.m));
end

function [ik,jk,k,n,m] = delayedkron(A,B)

[ma,na] = size(A);
[mb,nb] = size(B);

[ia,ja,sa] = find(A);
[ib,jb,sb] = find(B);
ia = ia(:); ja = ja(:); sa = sa(:);
ib = ib(:); jb = jb(:); sb = sb(:);
ka = ones(size(sa));
kb = ones(size(sb));
t = mb*(ia-1)';
ik = t(kb,:)+ib(:,ka);
t = nb*(ja-1)';
jk = t(kb,:)+jb(:,ka);

k=sb*sa.';
n=ma*mb;
m=na*nb;