%% function [util,f,w,prices,e] = TDMAweights(B,rates)%% Arguments:% B - routing matrix% rates - rates of links in different slots%% Output:% util   - optimal utility% f      - rates of flows% w      - weigths (durations) of slots% prices - shadow prices of links: a link with a high shadow price should be scheduled more often% e      - e != 0 if error (see fmincon).%% Description:%% find weights that maximizes utility fairness% we have the following system, when linearizing utility u_i:%% max 	sum u_i% subj. 	B*f <= x% 			x = R * w%        u_i <= a_ij*f_i + b_ij   (for all i,j)%% that is B*f - R*w <= 0, or [B|-R]*[f/x] <= 0   (1)% %function [util,f,w,prices,e] = TDMAweights(B,rates)% extract vector dimensionsl = size(rates,1);          % number of linksm = size(rates,2);          % number of slotsnf = size(B,2);             % number of flows% linearize utilitymaxR = max(reshape(rates,l*m,1));positiveR = reshape(rates,l*m,1);positiveR = positiveR(find(positiveR > 0));minR = min(positiveR);% sometimes maxR is very high and leads to numerical errors if using it.% we use mean insteadmeanR = mean(reshape(rates,l*m,1));% HEURISTIC: if maxR and minR are too close, results are bad.% hence, if necessery, we linearize on a larger regionif maxR/minR < 1e10    [Au, bu, umin] = util([maxR*sqrt(1e10*minR/maxR),minR/sqrt(1e10*minR/maxR)], 10);else    [Au, bu, umin] = util([maxR,minR], 10);end% I gave up using meanR since it is sometimes very inaccurate!%[Au, bu, umin] = util(meanR, 10);% [u,f,alpha]x0 = [zeros(nf,1)+umin;zeros(nf,1);ones(m,1)/m];Ae = [zeros(1,2*nf),ones(1,m)];be = 1;A = [zeros(l,nf), B, -rates];               % ineq. (1)A = [A;zeros(m,2*nf),diag(-ones(m,1))];b = zeros(l+m,1);for j=1:size(Au,1)	A = [A; diag(ones(nf,1)), -Au(j)*diag(ones(nf,1)), zeros(nf,m)];	b = [b; bu(j)*ones(nf,1)];endlb = [-Inf * ones(nf,1); zeros(m+nf,1)];ub = [];ff = [-ones(nf,1); zeros(nf+m,1)];options = optimset;options.LargeScale = 'off';options.MaxIter = 1000;options.TolFun  = 1.00e-06;options.Display = 'off';[x,util,e,output,lambda] = linprog(ff,A,b,Ae,be,lb,ub,x0,options);util = -util;u = x(1:nf);f = x(nf+1:2*nf);w = x(2*nf+1:2*nf+m);prices = lambda.ineqlin(1:l);