function [H,loss]=soc_avg(Gy,Y)
% 
if (nargin == 0) 
    disp('usage  [H,loss]=soc_avg(Gy,Y)')
   return
end

[ny,nu]=size(Gy);
 Jn = []; Fn = []; GnyT = [];
 Juu = eye(nu,nu);
 Juuhalf = sqrtm(Juu);
    for i = 1:nu
        GnyT = daug(GnyT,Gy');
        Jn = [Jn;Juuhalf(i,:)'];
        Fn = daug(Fn,Y);
    end
    Fxb = Fn*Fn';
    F = 2*Fxb;
    c = zeros(size(F,1),1);
    options = optimset('quadprog');
    options = optimset('TolFun',1e-10);    
    options = optimset('MaxIter',1000);
    [xn, fVal, ExitFlag, Output, lam] = quadprog(F,c,[],[],GnyT,Jn,[],[],[],options);
      H = reshape(xn,ny,nu)';
     loss = norm((Juuhalf*inv(H*Gy))*(H*Y),'fro');

