function [exitFlag exitMsg netStruct] = fa_minNetDelay_xde (netStruct, algorithmParameters, net2planParameters)

	% Algorithm for flow assignment problems which tries to minimize
	% the average network delay with the approximation for Poisson traffic. It uses a
	% flow-link formulation implemented in CVX. Av delay is given by:
	% (1/sum_d h_d) \sum_e y_e T_e
	% The delay in a link T_e is given by:
	% T_e = L_bits / (u_e - y_e) + t_e
	% Where:
	% h_d: traffic offered by demand d in bps
	% y_e: traffic traversing link e in bps
	% u_e: capacity of link e in bps
	% t_e: propagation delay of link e in seconds (net2planParameters needed to compute this)
	% L_bits: average packet length in bits (from net2planParameters)
	% 
	% It requires CVX solver installed and running.

	% Sanity check
	assert(isfield(netStruct, 'nodeXYPositionTable'), 'A set of nodes is required');
	assert(isfield(netStruct, 'linkTable'), 'A set of links is required');
	assert(isfield(netStruct, 'demandTable'), 'A set of demands is required');
	assert(isCvxInstalledAndRunning(), 'CVX solver must be installed and running');
		
	% initialize some variables
	L_bits = net2planParameters.averagePacketLengthInBytes * 8;
	R = net2planParameters.binaryRateInBitsPerSecondPerErlang;
	u_e = netStruct.linkCapacityInErlangs;
	t_e = netStruct.linkLengthInKm / net2planParameters.propagationSpeedInKmPerSecond;
	h_d = netStruct.offeredTrafficInErlangs;

	N = size(netStruct.nodeXYPositionTable, 1);
	E = size(netStruct.linkTable, 1);
	D = size(netStruct.demandTable, 1);
	
	% compute the adjacency matrix and the right-hand-side of flow
	% conservation constraints
	A_ne = linkTable2incidenceMatrix(netStruct.linkTable, N);
	w_nd = linkTable2incidenceMatrix(netStruct.demandTable, N);
	   
	cvx_begin

		cvx_quiet(true);
		variables rho_e(1,E); % utilization in link e 
		variables x_de(D, E); % fraction [0...1] of traffic of d, that traverses link e

		% NOTE: It is better to use rho_e instead of y_e as variable, since
		% rho_e ranges 0...1, and this helps to limit numerical errors
		
		% first sum is the contribution of queueing and transmission delays
		% second sum is the contribution of propagation delay
		minimize( (1 / sum(h_d)) * ( (L_bits/R*sum(quad_over_lin(rho_e, rho_e - rho_e.^2, 0))) + sum(rho_e .* u_e .* t_e) ) );

		subject to

			% flow conservation constraints
			% sum(e out n) x_de - sum(e in n) x_de == 1 if n=a(d)
			% sum(e out n) x_de - sum(e in n) x_de == -1 if n=b(d)
			% sum(e out n) x_de - sum(e in n) x_de == 0 otherwise
			% for all demand d, node n
			% a. vectorial form (all constraints in a single line of code)
			A_ne * x_de' == w_nd;
			% b. see fa_minimaxUtilization_xde for an equivalent expression with for-loops (slower)
		
			% define rho_e as te utilization of the link
			% rho_e = sum(d) h_d x_de / u_e, for all e
			rho_e == h_d * x_de ./ u_e; 
  
			% traffic in a link is <= its capacity (rho_e <= 1)
			rho_e <= 1;

			% non negativity
			x_de >= 0;
			% x_de <= 1; % this is redundant, so can be removed
	cvx_end

	if strcmp(cvx_status, 'Solved') || strcmp(cvx_status, 'Inaccurate/Solved')
			
		exitFlag = 0;
		exitMsg = 'OK!';
			
		% limit numerical errors (somthing close 0 is 0, close to 1 is 1)
		x_de(x_de > 1-net2planParameters.precisionFactor) = 1;
		x_de(x_de < net2planParameters.precisionFactor) = 0;

		% convert the x_de vector to a x_dp routing
		[x_dp pathList] = xde2xdp(x_de, netStruct.linkTable, ones(1, E), netStruct.demandTable);

		netStruct.routingMatrix = x_dp;
		netStruct.pathList = pathList;

	else

		exitFlag = -1;
		exitMsg = cvx_status;
		netStruct = [];

	end
		

end