function [] = demo() clc; clear all; %% Market share draft Lambdaiwant = 1.501; % Number of points npoints = 100; % Number of muN settings nmuN = 1; MSfast = zeros(nmuN,npoints); MSfastub = zeros(nmuN, npoints); MSfastlb = zeros(nmuN, npoints); Lambdavalues = zeros(nmuN, npoints); muNsettings = [1.1]; settingnumber = 0; for muN=muNsettings settingnumber = settingnumber+1 a = 0.5; low = 0.05; high = (1 + muN)*0.94; stepsize = (high - low)/(npoints - 1); count = 0; % Original Tmin Tmin = 10; for Lambda = low : stepsize : high count = count + 1 [pis, iMax, jMax, Tmin] = StationaryProbabilities( muN, Lambda, a, 20); [lambdaslow, lambdafast] = FindArrivalRates( Lambda, iMax, jMax, muN, a, pis ); Lambdavalues(settingnumber,count) = Lambda; MSfast(settingnumber,count) = lambdafast/Lambda; end; end; end function [ pis, iMax, jMax, T] = StationaryProbabilities( muN, Lambda, a, Tmin) %% This script finds the stationary probabilities when both A and N announce. Without loss of generality, let \mu_A = 1 and \mu_N take a higher value % T is the threshold that we move to attain an epsilon condition % pis is a matrix of stationary probabilities imax X jmax. i indexes slow server A, j % indexes fast server N. % Direction = -1 implies muN is barely less than its value. +1 implies muN is % barely more than its value. if(Lambda/(1+muN) > 0.95) epsilon = 1e-3; else epsilon = 1e-4 ; end; % What is the minimum I should allow T to be? % Tmin = 20; %Lambda = Lambda - 0.0001; % Initialize condition = 1; T = Tmin-1; % This would go into a loop while(condition > epsilon) T = T+1; iMax = T; jMax = ceil(T * muN); % Solve Ax = b A = zeros((iMax+1) * (jMax+1), (iMax+1) * (jMax+1)); b = zeros((iMax+1) * (jMax+1),1); % Put in the coefficients for A - one for each state (except 0,0) and one % for normalization eqcount = 0; for i=0:iMax for j=0:jMax % Ignore state (0,0) if(i==0 && j==0) continue; end; % Otherwise eqcount = eqcount + 1; b(eqcount) = 0; % Find the coefficients, which depend on the values of i and % j if(i==0) A(eqcount,Scalarize(i,j,jMax+1)) = Lambda + muN; A(eqcount,Scalarize(i+1,j,jMax+1)) = -1; if(j==1) A(eqcount,Scalarize(0,0,jMax+1)) = -Lambda * a; end; if(j<muN*T) A(eqcount,Scalarize(i,j+1,jMax+1)) = -muN; end; continue; end; if(j==0) A(eqcount,Scalarize(i,j,jMax+1)) = Lambda + 1; A(eqcount,Scalarize(i,j+1,jMax+1)) = -muN; if(i==1) A(eqcount,Scalarize(0,0,jMax+1)) = -Lambda *(1-a); end; if(i<T) A(eqcount,Scalarize(i+1,j,jMax+1)) = -1; end; continue; end; if(i==T && j==muN*T) A(eqcount,Scalarize(i,j,jMax+1)) = 1 + muN; A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda; continue; end; if(j==muN*T) A(eqcount,Scalarize(i,j,jMax+1)) = Lambda + muN + 1; A(eqcount,Scalarize(i+1,j,jMax+1)) = -1; A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda; if(j-1 == muN * i) A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda*a; elseif(j-1 < muN * i) A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; end; %{ if(direction == -1) if(j-1 <= muN*i) A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; end; elseif(direction == 1) if(j-1 < muN*i) A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; end; end; %} continue; end; if(i==T) A(eqcount,Scalarize(i,j,jMax+1)) = Lambda + muN + 1; A(eqcount,Scalarize(i,j+1,jMax+1)) = -muN; A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; if(j == muN * (i-1)) A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda*(1-a); elseif(j > muN * (i-1)) A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda; end; %{ if(direction == -1) if(j >= muN*(i-1)) A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda; end; elseif(direction == 1) if(j > muN*(i-1)) A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda; end; end; %} continue; end; % Otherwise, it is a generic state (i,j) A(eqcount,Scalarize(i,j,jMax+1)) = Lambda + muN + 1; A(eqcount,Scalarize(i+1,j,jMax+1)) = -1; A(eqcount,Scalarize(i,j+1,jMax+1)) = -muN; if(j > muN*(i-1)) A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda; elseif (j == muN*(i-1)) A(eqcount,Scalarize(i-1,j,jMax+1)) = -(1-a) * Lambda; end; if(j-1 < muN*i) A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; elseif(j-1 == muN * i) A(eqcount,Scalarize(i,j-1,jMax+1)) = -a * Lambda; end; %{ if(direction==-1) if(j >= muN*(i-1)) A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda; end; if(j-1 < muN*i) A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; end; elseif(direction==1) if(j > muN*(i-1)) A(eqcount,Scalarize(i-1,j,jMax+1)) = -Lambda; end; if(j-1 <= muN*i) A(eqcount,Scalarize(i,j-1,jMax+1)) = -Lambda; end; end; %} end; end; % Add the normalization equation eqcount = eqcount + 1; A(eqcount,:) = 1; b(eqcount) = 1; x = A\b; % Squarify pis = vec2mat(x,jMax+1); condition = sum(pis(iMax+1,:))+sum(pis(:,jMax+1)); %condition = pis(iMax+1, jMax+1); end; %T end
We use cookies to provide and improve our services. By using our site, you consent to our Cookies Policy. Accept Learn more