SDPSlover.rar

% SDP Solver clc clear tic; %SDPA格式 data=[ 2 ... 2 ... 2 2 ... 10.0 20.0... 0 1 1 1 1.0... 0 1 2 2 2.0... 0 2 1 1 3.0... 0 2 2 2 4.0... 1 1 1 1 1.0... 1 1 2 2 1.0... 2 1 2 2 1.0... 2 2 1 1 5.0... 2 2 1 2 2.0... 2 2 2 2 6.0]; data = data'; global n ns nis n2 n2s n2is nt nts ntis nblocks; % m = 约束矩阵数量 m = data(1); % nblocks = 约束矩阵的块对角线结构中的块数 nblocks = data(2); % 1D 索引 % n's = 每个nxn块的1D大小(n)(0是前缀) ns = data(3:3+nblocks-1); ns = [ 0; abs(ns)]; % n = sum( n's ) n = sum(ns); % n_i's = 每个块的累积一维指数 nis = zeros(size(ns)); for i = 1:length(ns) nis(i) = sum(ns(1:i)); end % 2D 索引 % n^2's = 每个nxn块的2D大小(n*n)(0前缀) n2s = ns.^2; % n2 = sum( n^2's ) n2 = sum(n2s); % n^2_i's = 每个块的累积2D索引 n2is = zeros(size(ns)); for i = 1:length(ns) n2is(i) = sum(n2s(1:i)); end % svec 索引 % nt's = (n*(n+1)/2) svec的三角索引(nxn块) (0前缀) nts = ns.*(ns+1)/2; % nt_i's = svec的累积三角索引(nxn块) ntis = zeros(size(nts)); for i = 1:length(nts) ntis(i) = sum(nts(1:i)); end % nt = sum(nt's) = sum( n_i*(n_i+1)/2 ) nt = sum(nts); % 从文件中读取数据并分配空间 b = -data(3+nblocks:3+nblocks+m-1); data = data(3+nblocks+m:end); c = zeros(n2, 1); A = zeros(m, n2); % 读取约束矩阵 for idx = 1:5:size(data,1)-1 matno = data(idx); %第几个约束矩阵 blkno = data(idx+1); %第几个块 i = data(idx+2); %行 j = data(idx+3); %列 entry = data(idx+4);%值 id1 = (j-1)*ns(blkno+1) + i + n2is(blkno); %行索引 id2 = (i-1)*ns(blkno+1) + j + n2is(blkno); %列索引 if matno == 0 c(id1) = -entry; c(id2) = -entry; else A(matno, id1) = -entry; A(matno, id2) = -entry; end end % 将A和c转换为svec版本 c = mat(c); c = svec(c); A_temp = A; A = zeros(m, nt); for i = 1:m A(i,:) = svec(mat(A_temp(i,:))); end % 初始X，Z选取，基于SDPT3文档 % 设ksi ksi = zeros(nblocks,1); temp = zeros(m, 1); for j = 1:nblocks for k = 1:m %求svec A中某个块的所有元素的范数，再用1+b/1+normA块， 获得缩放系数ksi temp(k) = (1 + abs(b(k))) / (1 + norm(A(k, vblock(j)))); end ksi(j) = ns(j+1) * max(temp); end % 设缩放系数eta eta = zeros(nblocks, 1); temp = zeros(m,1); for j = 1:nblocks for k = 1:m temp(k) = norm(A(k, vblock(j)));%求svec A中某个块的所有元素的范数 end eta(j) = (1/sqrt(ns(j+1))) * (1 + max( max(temp), norm(c(vblock(j))))); end % 设置初始迭代使用的X和Z X = zeros(n); for j = 1:nblocks X(block(j)) = eye(ns(j+1))*ksi(j); end x = svec(X); y = zeros(m,1); Z = zeros(n); for j = 1:nblocks Z(block(j)) = eye(ns(j+1))*eta(j); end z = svec(Z); % 计算每个块的normA，保存每个块的sqrt normsA = zeros(nblocks,1); for j = 1:nblocks normsA(j) = norm(A(:, vblock(j)), 'fro');%矩阵的各个元素绝对值的最大值 end normsA = max(1, sqrt(normsA)); % 计算normb normb = max(1, norm(b)); % 计算每个块 normc , 保留max normc = 1; for j = 1:nblocks normc = max(normc, norm(c(vblock(j)))); end % 缩放约束矩阵与系数向量 for j = 1:nblocks A(:, vblock(j)) = A(:, vblock(j)) / normsA(j); c(vblock(j)) = c(vblock(j)) / (normc*normsA(j)); end b = b / normb; % 缩放初始迭代x和z for j = 1:nblocks x(vblock(j)) = x(vblock(j)) * normsA(j); z(vblock(j)) = z(vblock(j)) / (normc*normsA(j)); end % 设置初始值 rp = b - A*x; Rd = c - z -A'*y; relgap = dot(x,z) / (1 + max(abs(dot(c,x)), abs(dot(b,y)))); %停止条件1：相对对偶间隙relgap pinfeas = norm(rp) / (1 + norm(b)); dinfeas = norm(Rd) / (1 + norm(c)); phi = max(pinfeas, dinfeas); %停止条件2：不可行系数phi soln_relgap = 1e-6; %停止准则 soln_phi = 1e-6; %停止准则 gamma = 0.9; iter = 0; % 更新中间解直到收敛 while(relgap > soln_relgap || phi > soln_phi) % 得到X和Z的Cholesky分解 X = smat(x); Q = chol(X); Z = smat(z); P = chol(Z); % 计算残差mu，relgap，phi %(4) rp = b - A*x; Rd = c - z -A'*y; %svec Rd mu = dot(x,z) / n; %mu relgap = dot(x,z) / (1 + max(abs(dot(c,x)), abs(dot(b,y)))); phi = max(norm(rp)/(1+norm(b)), norm(Rd)/(1+norm(c))); Rc = -svec( H(X*Z, P) ); % 计算预测方向（伪方向 sigma =0） (10-16) HKM方向 得到dx，dz的svec形式 M = A*skmult(X, Z\eye(n), A'); %NT （43） h = rp + A*(skmult(X, Z\eye(n), Rd) - skmult(P\eye(n), P\eye(n), Rc)); dy = M\h; dz = Rd - A'*dy; %dx = -x - skmult(X, Z\eye(n), dz); dx = skmult(P\eye(n), P\eye(n), Rc) - skmult(X, Z\eye(n), dz); % 计算 alpha (原始步长)(33) alphas = zeros(1,nblocks); for j = 1:nblocks XDX = (Q(block(j))') \ smat_1block(dx(vblock(j))) / Q(block(j));%（32） eigs_a = sort(eig(XDX)); %排序 lambda_a = eigs_a(1); %取最小的特征值 if lambda_a < 0 alphas(j) = -1/lambda_a; else alphas(j) = Inf; end end alpha = min(1, gamma*min(alphas)); % 计算 beta (对偶步长)(34) betas = zeros(1,nblocks); for j = 1:nblocks ZDZ = (P(block(j))') \ smat_1block(dz(vblock(j))) / P(block(j)); eigs_b = sort(eig(ZDZ)); lambda_b = eigs_b(1); if lambda_b < 0 betas(j) = -1/lambda_b; else betas(j) = Inf; end end beta = min(1, gamma*min(betas)); % 计算指数e (sigma的指数) if mu > 1e-6 if min(alpha, beta) < 1 / sqrt(3) e = 1; else e = max(1, 3*min(alpha, beta)^2); end else e = 1; end % 计算 sigma (中心参数) if dot(x+alpha*dx, z+beta*dz) < 0 sigma = 0.8; else frac = dot(x+alpha*dx, z+beta*dz)/dot(x,z); sigma = min(1, frac^e); end % 计算 正确 的预测方向 (10-16) HKM方向 得到dx，dz的svec形式 Rc = svec(sigma*mu*eye(n)-H(X*Z, P)-H(smat(dx)*smat(dz), P)); h = rp + A*(skmult(X, Z\eye(n), Rd) - skmult(P\eye(n), P\eye(n), Rc)); dy = M\h; dz = Rd - A'*dy; dx = skmult(P\eye(n), P\eye(n), Rc) - skmult(X, Z\eye(n), dz); gamma = 0.9 + 0.09*min(alpha, beta); % 计算正确 alpha (原始步长)(33) for j = 1:nblocks XDX = (Q(block(j))') \ smat_1block(dx(vblock(j))) / Q(block(j)); eigs_a = sort(eig(XDX)); lambda_a = eigs_a(1); if lambda_a < 0 alphas(j) = -1/lambda_a; else alphas(j) = Inf; end end alpha = min(1, gamma*min(alphas)); % 计算正确 beta (对偶步长)(34) for j = 1:nblocks ZDZ = (P(block(j))') \ smat_1block(dz(vblock(j))) / P(block(j)); eigs_b = sort(eig(ZDZ)); lambda_b = eigs_b(1); if lambda_b < 0 betas(j) = -1/lambda_b; else betas(j) = Inf; end end beta = min(1, gamma*min(betas)); % 更新步骤 x = x + alpha*dx; y = y + beta*dy; z = z + beta*dz; gamma = 0.9 + 0.09*min(alpha, beta); % 打印迭代结果 iter = iter + 1; primal = dot(c,x); %原始结果 dual = dot(b,y); %对偶结果 pinfeas = norm(rp) / (1 + norm(b)); %原始可行性 dinfeas = norm(Rd) / (1 + norm(c)); %对偶可行性 fprintf('%d\t%.3f\t%.3f\t%.3e\t%.3e %.3f\t%.2e %.2e %.2e %.2e %.2e %.2e %.2f %2e %2e\n', ... iter, alpha, beta, sigma, mu, e, pinfeas, dinfeas, relgap, primal, dual, phi, gamma, norm(rp), norm(Rd)); end disp(' '); disp('计算完毕'); % 反缩放约束矩阵 for j = 1:nblocks A(:, vblock(j)) = A(:, vblock(j)) * normsA(j); c(vblock(j)) = c(vblock(j)) * (normc*normsA(j)); end b = b * normb; % 反缩放结果 for j = 1:nblocks x(vblock(j)) = x(vblock(j)) * normb / normsA(j); z(vblock(j)) = z(vblock(j)) * (normc*normsA(j)); end y = y * normc; X = smat(x); Z = smat(z); % 打印结果 disp(['原始结果: ', num2str(dot(-c,x))]); disp(['对偶结果: ', num2str(dot(-b,y))]); time = toc; disp(['时间: ', num2str(time), ' s']);