The most important contribution of this study is that we extends the original work of L
& Z. We fit it into an SC domain by the adequate “transformation” of non-tandem systems.
The sophisticated processing logic inside each echelon can also be handled by our proposed approach, for example server unavailability herein. Though not exhibited herein, we believe this approach can also handle any arbitrary service distribution since phase type can approximately represent any probability distribution (Svoronos and Zipkin, 1991) and
⋅/PH/1, where the arrival process can be of Markovian or phase-type, is a well-developed domain of QBD process (Neuts, 1994).
We also show the possibility of extending the original model to handle other order/replenish scheme. Specifically that of (r, q) type is successfully tested though the accuracy is satisfactory only under specific conditions. Integrating the results of QBD modeling for parallel processing with L & Z is also tested with limited success. To make L
& Z more general, more efforts are needed in the future.
In conclusion, we have demonstrated that by using the matrix analytic approach, the evaluation of a complex SC using base stock as control logic, performs as expected through simulation verification. The relative errors between Approximation and Simulation are all below 10% for retailers adopting MTO policies in our test problems. When all the retailers also adopt the MTS policy, numerical studies show that the approximation is accurate for medium traffic intensity and acceptable for high traffic intensity. This shows that the matrix analytic approach by combining L & Z and other decomposition technique such as QBD for supply chain analysis is feasible, therefore its application is not limited to tandem processing queue as reported by L & Z and Zipkin (1995) but also for broader context. Herein we show its application on a tandem SC where the end stage is a distribution system with identical or non-identical retailers.
102
The non-identical case refers to non-identical demand rates with non-identical distribution processing rates. Also we transformed and incorporate a fork-join type supplier’s process into our evaluation model to make it more general. In the literature of the stochastic production-distribution system, most models are developed and analyzed separately. Those models are usually difficult to integrate as one single model.
Next, the MMPP models are developed for different scenarios related to uncertainty to investigate respective system behavior for each individual stage under uncertainty concern.
After linking together all the stages the system behavior of the whole chain is analyzed by queueing network analysis (QNA) under MTO supply mode. However, the QBD decomposition for parallel processing under various uncertainties may not be suitable under MTS mode. From the study herein, the application only works under specific conditions. The application of the analytic approach on resource allocation, specifically server and repairman decision is demonstrated through an illustrative example. The results show the analysis provides valuable managerial insights.
Finally, we demonstrated various optimization methods for the optimization of various SC topologies under service constraints. We showed the impact on system performance due to upstream unavailability or poor quality under different service levels. The unreliability concern in SCM is seldom explored in the literature, to our knowledge.
Simchi-Levi et al. (2003) mentioned that the impact on the performance of an SC due to 911-terror-attack in 2001 is influential, including delayed delivery, delayed customs, poor communications. On the other hand, the impact on an SC due to supply and demand uncertainty in the long run is obscure. Herein, we assessed the above argument through quantitative analysis. To sum up, we characterize the major contributions of this work:
1. Extending L & Z to allow not just for tandem supply network but also for more general supply form. The impact of upstream unavailability on the performance of an SC in the long run is made clear through empirical study.
103
2. Proposing feasible modeling approach to investigate the impact of non-stationary demand and unreliable service processes on the performance of an SC.
3. Relaxing the single server assumption usually used in tandem queueing models for SC analyses. Though numeric study is satisfactory for low traffic intensity only.
4. Exploring the possibility of extending the analytic model to allow for other inventory control policies.
5. Investigating cases in applying classical or modern optimization methods for solving stochastic SC problems. The obtained results consolidate several previous researches with new insights.
We also prove that the average number of operative machines is equal (proportional) to the average number of machines under repair when mean time to failure and mean time to repair are the same (proportional) by using a matrix algebraic approach. This property is not explicitly related to repairman. During the process of this study we encountered several problems, which raised opportunities for further study. Below we itemize these problems and propose some possible strategies for research.
Evaluation model
The QBD combined with L & Z forms the mainframe of our evaluation model. We tested the performance of the proposed evaluation model in a tandem setting. As can be seen herein, the approximation model has its limitation, to attain a more general setting, alternative evaluation model consisting of new decomposition method and new approximation method may be necessary for other topology and system dynamic of an SC not addressed herein. For example, problems dealing with generalized QN with arbitrary input and service rate and planned inventory setting are still open.
System dynamic inside each echelon
In this work we used abstract level of SC modeling. Though chapter 4 discussed several details in modeling the behavior of an SC, it’s not enough. To unveil the System
104
dynamic inside each echelon and relate it to system performance, more thorough studies should be done. For example, herein we assumed identical processing rates inside assembly logic, what if they are non-identical. Also we used simple model to represent transportation activity. However, practical logistics is usually more complicated than the model formed herein. Optimizing the SC model with more practical logistics concern is an interesting area. Also future study may investigate the impact of uncertainties on resource design decision as functions of correlated effect of MMPP input.
Optima seeking
Due to accuracy reason, we do not put all the parameters investigated herein into a nutshell and optimize it. However, such modeling approach may be necessary for practical use. Further, the evaluation function of the studied problem is more complicated then other engineering optimization problems using the same meta-heuristic search methods. When evaluation model is changed or becomes more intractable because of more complex SC structure and/or system dynamics, more efficient and effective search algorithm other than meta-heuristic as suggested herein may have to be developed. Heuristic approaches such as those in Boyaci and Gallego (2001) are examples. Alternatively, in addition to GA and SA, other state-of-the-art meta-heuristics techniques such as Ant colony algorithm, scatter search, Tabu search among others, may provide different solution flavors.
Control policy
As we reviewed in chapter 2, there are many other control policies which are attributed to “pull type” production/inventory control as studied herein. In this study we used this (r, q) control policy by analogous transformation of the performance formula of Zipkin (2000).
It works well under some situations. However, further study regarding extension of the analytic model to other inventory control schemes, such as (s, S) or KANBAN may have to be addressed in the future. On the other hand, periodic review policy as opposed to continuous review policy studied herein may also be investigated.
105
Buffer design
We used dual-buffer design herein to develop our evaluation model. The input buffer is assumer infinite. There are more SC designs with only one finite buffer due to space limitation or cost concern. The single buffer design has to take blocking effect into its modeling logic. It may be interesting to compare the performance between single and dual buffer designs under various SC topologies. Also dual buffer with finite waiting line is an alternative design, which may be more realistic due to practical reason.
Value of the information
In this study, we used centralized information sharing scheme. When demand is generated at the end stage, all of the stages along the chain immediately reflect this information. It is well known that decentralized information sharing scheme is the cause of so-called “bullwhip effect” in an SC. The demand information is exaggerated toward upstream stages and the variation of demand is amplified along the chain to upstream stages. Further study may investigate the saving of cost due to the centralized information.
Push-pull boundary
In Simchi-Levi et al. (2003), “pull” type of inventory control policy is suitable for the case when demand is highly uncertain and no economy of scale exists. Under these conditions MTO supply mode is adequate. On the contrary, when demand uncertainty is low and economy of scale exists, MTS supply mode is more adequate. However, mixed type of supply mode is possible when one adequately selects the so-called “Push-pull boundary”. Before the boundary, push (MTS) supply is conducted while pull (MTO) supply is executed after the boundary. The performance comparison between our model and the design of this mixed policy is an interesting research area.
The role of the warehouse
Recently, there are several arguments regarding the role of the warehouse in distribution stage. Cross-docking policy argues the inventory holding cost can be largely
106
saved and so is the supply lead-time. Direct sale policy argues there is no need for the existence of the warehouse. The supply of the manufacturer may directly serve the end stage retailers to save operation cost. Since it is obviously shown herein, the complex interaction of the SC decides the performance of an SC, more thorough studies may have to be conducted to justify which policy is most suitable under what conditions.
As we stated in literature review of chapter 2, high-level modeling approach may be another direction for SCM such as Stochastic Petri net (SPN). Most researches on SPN or generalized SPN (GSPN) can be found in computer and telecommunication areas.
Alexander (2001) first showed how to use MGM as the main solving process for QBD model. They illustrated how to map from an infinite Stochastic Petri net (iSPN) model to a QBD process for MGM analysis. Trivedi (2002) illustrated the technique of transforming the reachability graph (RG) obtained from SPN to corresponding CTMC generator matrix in detail. Arns et al. (2002) used Proc/B, a notational description tool, which can then be translated to queueing network or SPN for performance calculation. They used the notation model to compare ordinary distribution and web retailing by constructing synchronized (parallel) sub-models of Proc/B as SPN model and solve the underlying CTMC and then aggregated them as a single queue. Then the overall QN is solved by product-form solution.
Under their approach, all the underlying transformation and calculation were automatically accomplished.
107
Appendix A QBD process and simulation model
A.1 QBD Process
Define the phase type distribution of a Markov process as a stochastic process having parameters m ,G* ,andα* if it can be expressed as a first-passage time random variable, then be represented in the following matrix form:
⎥⎥
Where Λ=λI . Then we can find the steady-state probability vector )
Combining (A.1.1) with the equations PQ=0 yields the following system of equations,
The solution of (A.1.2) involves the characteristic equation
* + =0.
+R(G -Λ) R2G∆α*
Λ (A.1.3)
108
The matrix-geometric solution of (A.1.3) is
,
whose dimension is chosen to fit the context.
A.2 Simulation Model of 3.4 by using Arena (Kelton et al., 2002)
Cust om er Ar r ive unit dem and war ehouse on hand
backor der
r eplenishm entwar ehouse Separ at e 3
m yopic allocat ion
Sta t i o n 1
Fill B2 Backor der
B2 Avail? Fill Cust om er
Backor der
109
Appendix B Moment derivations and related proofs
B.1 First and second moment derivations of chapter 4
Here we show the derivation of u(1), which is not seen in Neuts and Lucanton (1979).
B.2 Proof of proposition 1
Assume the generator of a particular server queue is as follows. All the sub-matrix is similar to (4.4) of chapter 4 with m servers.
110
Since πQ = , after eliminating the last equivalence equation, we have 0 0 Substitute (B.2.3) into (B.2.2) and after eliminating and arrange terms we have
2
Continue doing in this way, finally we get
1 .
Multiply by unity vector of proper dimension on both sides yield
1
.
Plug in the relationship
111
.
and ) 1
( ,
i i i
J i γ
ζ
=
=
Γ
−
Ψ I
If γ = , from the definitions of E[O] and E[RE], we have the result. ζ B.3 Proof of corollary 1
The proof immediately follows from PROPOSITION 1 by rearranging terms.
112
Appendix C Matlab codes of algorithms
C.1 CSA
CSA contains 3 functions: CSABCAP, simannBcap, OptSABCap.
C.1.1 CSABCAP code function CSABCAP global N lamda h1 h2 fid fid = fopen('CSAout.m','w');
% Initialize t=cputime;
sa_t=85;
sa_rt=0.8;
sa_nt=5;
sa_ns=20;
fun_name='OptSABCap';
run=5;
s_l=0.1;
for N=4:7
for iter=1:run
c=cell(1,N);p=c;r=c;pai=c;
B=zeros(1,N);W=B;I=B;h1=B;h2=B;lamda=1;delta=0.1;
i=1:N;h1=0.5+delta*(i-2);h2=0.5*(1+delta).^(i-1);
u=zeros(1,N);v=u;
for i=1:N
c{i}=zeros(i);p{i}=c{i};r{i}=zeros(1,i);pai{i}=r{i};
end
% Decide searching bounds
UB1=100*ones(1,N);LB1=zeros(1,N);
UB2=100*ones(1,N);LB2=1.01*ones(1,N);
% Run CSA
[xopt,fopt,SLopt,rhopt,Wopt,Iopt,Bopt]=SimAnnBcap(fun_name, LB1, UB1,…
LB2, UB2, sa_t, sa_rt, sa_nt, sa_ns, s_l);
fprintf(fid,'N %d Run %d xopt, fopt, SLopt, rhopt, Wopt, Iopt, Bopt\n',N,iter);
fprintf(fid,'(');
for i=1:N if i==N
fprintf(fid,'%d)',xopt(i));
else
fprintf(fid,'%d, ',xopt(i));
end end
for i=N+1:2*N
fprintf(fid,'%.2f ',xopt(i));
end
fopt_1=fopt-sum(xopt(N+1:2*N));
fprintf(fid,' %.3f %.3f %.3f ',fopt,fopt_1,SLopt);
fprintf(fid,'(');
for i=1:N
113
if i==N
fprintf(fid,'%.2f)',rhopt(i));
else
fprintf(fid,'%.2f, ',rhopt(i));
end end for i=2:N
fprintf(fid,'%f ',Wopt(i));
end for i=1:N
fprintf(fid,'%f ',Iopt(i));
end for i=1:N
fprintf(fid,'%f ',Bopt(i));
end
fprintf(fid,'\n');
end end
fclose(fid);
e=cputime-t;
fprintf('operation duration: %f',e);
C.1.2 simannBcap code
function [xopt,fopt,SLopt,rhopt,Wopt,Iopt,Bopt]=simannBcap(func, LB1, UB1, LB2, … UB2, sa_t, sa_rt, sa_nt, sa_ns, s_l)
global N R fid
sa_neps=4; %number of times eps sa_eps=1e-6; %convergence criteria
sa_maxeval=500000; %maximum number of function evaluations sa_nacc=0; %number of acceptions
sa_nevals=0; %number of evaluations fstar=Inf*ones(sa_neps,1); %last optimum at each sa_nt x(1:N-1+R)=LB1(1:N-1+R)+(UB1(1:N-1+R)-LB1(1:N-1+R)).*rand(1, N-1+R);
%starting values for buffer
u(1:N-1+R)=LB2(1:N-1+R)+(UB2(1:N-1+R)-LB2(1:N-1+R)).*rand(1, N-1+R);
%starting values for capacity
[f,SL]=feval(func,round(x),u); %function evaluation with parameters x u xopt=[round(x) u];
fopt=f;
SLopt=SL;
fstar(1)=f;
VM1=(UB1-LB1); %maximum step size VM2=(UB2-LB2); %maximum step size while 1
nup=0; %number of uphill movements nrej=0; %number of rejections
nnew=0; %number of new global optimum ndown=0; %number of downhill movements lnobds=0;
114
nacp=zeros(N-1+R,1);
for m=1:sa_nt % loop number for each temperature VM(step adjusting) occur for each
% loop
for j=1:sa_ns % loop number for the purpose of adjusting step, each variable try sa_ns
% times
for h=1:N-1+R % for each variable if sa_nevals>=sa_maxeval
disp('max function evaluations achieved') return
end xp=x;
up=u;
xp(h)=x(h)+VM1(h)*(2*rand(1,1)-1.0);
if (xp(h)<LB1(h)) | (xp(h)>UB1(h))
xp(h)=LB1(h)+(UB1(h)-LB1(h))*rand(1,1);
end
up(h)=u(h)+VM2(h)*(2*rand(1,1)-1.0);
if (up(h)<LB2(h)) | (up(h)>UB2(h))
up(h)=LB2(h)+(UB2(h)-LB2(h))*rand(1,1);
end
[fp,SL,rho,W,I,B]=feval(func,round(xp), up);
sa_nevals=sa_nevals+1;
if (fp<=f & SL<=s_l)
x=xp; % trial replace best sol x
u=up; % trial replace best sol u f=fp;
sa_nacc=sa_nacc+1;
nacp(h)=nacp(h)+1;
ndown=ndown+1;
if fp<fopt
xopt=[round(xp) up];
fopt=fp;
SLopt=SL;
Wopt=W;
Iopt=I;
Bopt=B;
rhopt=rho;
sa_opteval=sa_nevals;
nnew=nnew+1;
end
else % function value increases
p=exp((f-fp)/sa_t);% Metropolis' criteria pp=rand(1,1);
if (pp<p & SL <=s_l) x=xp;
u=up;
f=fp;
sa_nacc=sa_nacc+1;
nacp(h)=nacp(h)+1;
nup=nup+1;
115
else
nrej=nrej+1;
end end end end
%adjust maximal step size vm sa_ns times have passed, each variable
% hascalculaterespective nacp(i) if ratio of accept of the variable (less function
% value) or prob. of escape is high (function value not so high) then enlarge the
% search area within the neighborhood of the variable; if ratio of accept is low then
% shrink search area within the neighborhood of the variable;
c=ones(N-1+R,1)*2;
for i=1:N-1+R
ratio=nacp(i)/sa_ns;
% Step length adjustment, Note if 0.4<=ratio<=0.6 no change happens if ratio>0.6
VM1(i)=VM1(i) * (1+c(i)*(ratio-0.6)/0.4);
VM2(i)=VM2(i) * (1+c(i)*(ratio-0.6)/0.4);
elseif ratio <0.4
VM1(i)=VM1(i) / (1+c(i)*((0.4-ratio)/0.4));
VM2(i)=VM2(i) / (1+c(i)*((0.4-ratio)/0.4));
end
if VM1(i)>(UB1(i)-LB1(i)) VM1(i)=UB1(i)-LB1(i);
end
if VM2(i)>(UB2(i)-LB2(i)) VM2(i)=UB2(i)-LB2(i);
end end
for i=1:N-1+R nacp(i) = 0;
end end
% check termination criteria for the current temperature, if the current optimal f (before
% changing temperature) less than global optimal (fopt) within eps, quit; notice fp (new
% function value) could become current optimal f because of Metropolis selection process
% fstar(1)=f; f, f(1)* is current optimum for this temperature quit = (((fstar(1)-fopt) <= sa_eps) & (SL<=s_l));
% guarantee current optimum is global optimum
% if within 4 times of temperature reduction current optimal f and any of previous
% temporary optimal difference within eps, quit. means dwindling if any(abs(fstar-f)>sa_eps)
quit=0;
end if quit
disp(['simulated annealing achieved termination after ', num2str(sa_nevals),' evals']);
return end
sa_t=sa_t*sa_rt;% reduce temperature fstar(2:4)=fstar(1:3);
% continue from current optimum
116
x=xopt(1:N-1+R);
u=xopt(N+R:2*(N-1+R));
f=fopt;
end %while
C.1.3 OptSABCap code
function [TC,SL,rho,W,I,B]=OptSABCap(x,u)
% Performance evaluation (Object function) for constrained SA (CSA) global N lamda h1 h2
N=4;delta=0.1;i=1:N;lamda=1;h1=0.5+delta*(i-2);h2=0.5*(1+delta).^(i-1);
v=u-lamda;
c{1}=[-v(1)];r{1}=[1];T=0;
for i=2:N
c{i}=[c{i-1} zeros(i-1,1)-sum(c{i-1},2);zeros(1,i-1) -v(i)];W(i)=lamda/v(i);
end for i=1:N
p{i}=lamda*inv(lamda*eye(i)-c{i});
end
pai{1}=r{1}*p{1};
for i=2:N
r{i}=[r{i-1}*p{i-1}^x(i-1) 1-r{i-1}*p{i-1}^x(i-1)*ones(i-1,1)];pai{i}=r{i}*p{i};
end for i=1:N
B(i)=pai{i}*p{i}^x(i)*inv(eye(i)-p{i})*ones(i,1);I(i) = x(i) - pai{i}*inv(eye(i) - p{i})* … ones(i,1) + B(i);
end
SL=r{N}*p{N}^x(N)*expm(c{N}*T)*ones(N,1);
WIP=sum(I+W)-I(N);
TC=sum(h1.*W+h2.*I)+10*B(N)+sum(u);
rho=1./u;
C.2 CGA
CGA contains 13 functions: CGA, initialize_r0, validate_r0, GeneticAlgorithm, doHomomorphMap, SimpleCrossover, ArithmeticCrossover, Xmutation, fineBoundaryResolution, delta, deltaInverse, sum_tSegments, Optga.
C.2.1 CGA code function CGA
global nG nPi nP nSC nAC nVar bounds ctr Pc Pm fid FINE_RESOLUTION … COARSE_RESOLUTION xmapped Xdummy r c p pai W h1 h2
% GA parameters run = 5;
nG = 100; % total number of generation nPi= 24; % population size (Fixed) nP = 12; % Mating pool size
nSC = 6; % number of simple crossover - each yields 2 children nAC = 6; % number of Arithmatic crossover - each yields 2 children Pc = 0.8; % probability of crossover
Pm = 0.01; % probability of mutation MAP_RESOLUTION = 5.0e-5;
117
FINE_RESOLUTION = 8;% feasible space searching COARSE_RESOLUTION = 15; %feasible space searching fid = fopen('CGAout.m','w');
for nVar=4:7
bounds=repmat([0 30],1,nVar);
% Initialize
c=cell(1,nVar);p=c;r=c;pai=c;t=cputime;
B=zeros(1,nVar);W=B;I=B;h1=B;h2=B;lamda=1;WIP=0;delta=0.1;
i=1:nVar;h1=0.5+delta*(i-2);h2=0.5*(1+delta).^(i-1);
u(i)=1.25;v(i)=u(i)-lamda;
for i=1:nVar
c{i}=zeros(i);p{i}=c{i};r{i}=zeros(1,i);pai{i}=r{i};
end
Xdummy=zeros(nPi,nVar);xmapped=Xdummy;
XA=zeros(nP,nVar);NM=zeros(nPi,nVar);
% Compute c, p of L & Z '92 c{1}=[-v(1)];
for i=2:nVar
c{i}=[c{i-1} zeros(i-1,1)-sum(c{i-1},2);zeros(1,i-1) -v(i)];W(i)=lamda/v(i);
end
for i=1:nVar
p{i}=lamda*inv(lamda*eye(i)-c{i});
end r{1}=[1];
pai{1}=r{1}*p{1};
for iter = 1:run% random runs
% Initial generated feasible reference point [ret,r0]=initialize_r0;
if (ret==1)
fprintf('initial feasible reference point : ');
for i=1:nVar
fprintf('%f ',r0(i));
end
fprintf('generated\n');
end
% initial random population x = -1 + rand(nPi,nVar) * 2;
% HomomorphMapping
[xmapped]=doHomomorphMap(x,r0);
[fval] = Optga(xmapped);
[Yvect,Ip]=sort(fval);
% Initial populated population
Xgen=zeros(nPi,nVar);Fgen=zeros(nPi);
for i = 1:nPi
Xgen(i,:)=x(Ip(i),:);
Fgen(i)=Yvect(i);
end
Xbest=Xgen(1,:);
Fbest=Fgen(1);
fprintf(fid,'Run %d ---GA initial---\n',iter);
fprintf(fid,'Fbest: %f\n',Fbest);
118
ctr = 0; % Initialize iteration counter while (ctr~=nG)
ctr = ctr + 1;
[x,f,Xgen,Fgen]=GeneticAlgorithm(@Optga,Xgen,Fgen,r0);
fval(ctr) = f;
% Preserving the best feasible solution if (fval(ctr) < Fbest)
Fbest=fval(ctr);
Xbest=round(x);
fprintf(fid,'Generation:%d, Current best value: %f\n',ctr,fval(ctr));
end end
[invalid,nlineq] = validate_r0(Xbest);
fprintf(fid,'N: %d Run: %d ---GA results---\n',nVar,iter);
fprintf(fid,'Xbest: ');
for i=1:nVar
fprintf(fid,'%d ',Xbest(i));
end
fprintf(fid,'\n');
fprintf(fid,'Fbest: %f\n',Fbest);
fprintf(fid,'SL: %f\n',nlineq+0.1);
end end
fclose(fid);
e=cputime-t;
fprintf(1,'operation duration: %f',e);
C.2.2 initialize_r0 code function [ret,r0]=initialize_r0 global nVar bounds
TRIES=10000;
for k = 0:99
for count = 1:10 * TRIES for i = 1:nVar
r0(i) = rand(1,1)*(bounds(2*i)- bounds(2*i-1)) + bounds(2*i-1);
end
[invalid] = validate_r0(r0);
if(~invalid)
fprintf('Valid reference point found!!\n');
fprintf('Starting CGA...\n\n');
ret=1;
return end end end
fprintf('Cannot find valid reference point.\n');
ret=0;
C.2.3 validate_r0 code
function [invalid,nlineq] = validate_r0(x) global r c p pai W nVar h1 h2
119
x=round(x);invalid = 0;
for i=2:nVar
r{i}=[r{i-1}*p{i-1}^x(i-1) 1-r{i-1}*p{i-1}^x(i-1)*ones(i-1,1)];pai{i}=r{i}*p{i};
end
for i=1:nVar
B(i)=pai{i}*p{i}^x(i)*inv(eye(i) - p{i})*ones(i,1); I(i) = x(i) - pai{i}*inv(eye(i) - … p{i})*ones(i,1) + B(i);
end i=1:nVar;
%SL=r{nVar}*p{nVar}^x(nVar)*expm(c{nVar}*0)*ones(nVar,1);
nlineq(1)=r{nVar}*p{nVar}^x(nVar)*expm(c{nVar}*0)*ones(nVar,1) - 0.1;
for i = 1:1
if( nlineq(i) > 0.0 ) invalid = 1;
return end
end
C.2.4 GeneticAlgorithm code
function [XCurbest,FCurbest,Xgen,Fgen]=GeneticAlgorithm(objfun,Xgen,Fgen,r0)
% select top 50% populations to propogate to the next generation,
% Offspring replace the top 50% populations for each crossover operator, then
% non-uniform mutation for whole population global nG nPi nP nSC nAC nDC nI nVar bounds fid
% Record best chromosome in the beginning of each generation;
XBest=Xgen(1,:);
Best=Fgen(1);
% Operator --- SimpleCrossover, ArithmeticCrossover.
[XSChild] = SimpleCrossover(Xgen);
[XAChild] = ArithmeticCrossover(Xgen);
i=1:2*nAC;
Xgen(i,:) = XAChild(i,:);
Xgen(2*nAC+i,:) = XSChild(i,:);
[XNMu]= XMutation(Xgen);
Xgen = XNMu;
[n m] = size(Xgen);
[xmapped]=doHomomorphMap(Xgen,r0);
[fval] = Optga(xmapped);
Xdummy=Xgen;
[Yvect,Ip]=sort(fval);
XCurbest=xmapped(Ip(1),:);
FCurbest=Yvect(1);
clear Xgen fval% clear memory for next generation Xgen=zeros(nP,nVar);
i = 1:nP;
% Elitist (if current best not superior to last generation, replace the last chromosome w/ the
% last best) if (FCurbest<Best)
Xgen(i,:)=Xdummy(Ip(i),:);
Fgen(i)=Yvect(i);
else
120
Xgen(i,:)=Xdummy(Ip(i),:);
Fgen(i)=Yvect(i);
Xgen(nP,:)=XBest;
Fgen(nP)=Best;
End
C.2.5 doHomomorphMap code
function [xmapped]=doHomomorphMap(x,r0)
global nVar nPi bounds MAP_RESOLUTION COARSE_RESOLUTION fid [n,m]=size(x);
infeasible = 0;
j=1:n;
ymax(j)=max(abs(x),[],2);
% find sVect y=zeros(nPi,nVar);
for i=1:nVar
y(:,i)=x(:,i)./ymax';
end i=1:nVar;
s = repmat(1/2*(bounds(2*i) + bounds(2*I - 1)),nPi,1) + repmat( 1/2*(bounds(2*i) - … bounds(2*i-1)),nPi,1).*y;
s_r0=s-repmat(r0,nPi,1);
for j=1:nPi
% Get tVect tCount = 1;
for i = 1:COARSE_RESOLUTION
tempVect = r0 + s_r0(j,:)*(i / COARSE_RESOLUTION);
[invalid,nlineq] = validate_r0(tempVect);
if(invalid ~= infeasible) infeasible = invalid;
tCount = tCount + 1;
end end
oddCount = mod(tCount,2);
lastIsInvalid = invalid;
[invalid,nlineq] = validate_r0(s(j,:));
boundaryIsInvalid = invalid;
% Populate regions
tVect(1) = 0.0;infeasible = 0;tCount = 2;
for i = 1:COARSE_RESOLUTION
tempVect = r0 + s_r0(j,:)*(i / COARSE_RESOLUTION);
[invalid,nlineq] = validate_r0(tempVect);
if(invalid ~= infeasible)
[tVect(tCount)] = fineBoundaryResolution(s_r0(j,:),i,r0,infeasible);
tCount = tCount + 1;
infeasible = invalid;
end end
if(oddCount & (~lastIsInvalid) & (~boundaryIsInvalid)) tVect(tCount) = 1.0;
elseif(oddCount & (~lastIsInvalid) & boundaryIsInvalid)
121
tVect(tCount) = fineBoundaryResolution(s_r0(j,:), COARSE_RESOLUTION, r0, … infeasible);
elseif((~oddCount) & lastIsInvalid & (~boundaryIsInvalid))
[tVect(tCount)]=fineBoundaryResolution(s_r0(j,:),COARSE_RESOLUTION,r0, ...
infeasible);
tCount = tCount + 1;
tVect(tCount) = 1.0;
else
tCount = tCount -1;
end
% Do mapping!!
[deltaInvVal] = deltaInverse(ymax(j), tVect, tCount);
xmapped(j,:) = r0 + s_r0(j,:)*deltaInvVal;
end
C.2.6 SimpleCrossover code
function [XS] = SimpleCrossover(X)
% creating children by simple 1-point crossover global nVar nSC Pc
first=0;one=0;
XS = X;
for i=1:2*nSC r = rand(1,1);
if (r<Pc)
if (r<Pc)