It seems that the only way to get rid of analytical solutions is to make the problem extremely hard such that nobody manages to derive the analytical solution, or the analytical solution simply does not exist...
The problem may be very simple. The computation cost of analytical solution (if it exist) has to be larger than the cost of monte carlo.
I agree, many Cody problems take the opposite route (they force players to find analytical solutions by using very stringent error tolerances or very large numbers/samples) but very few problems favor algorithmic approaches over analytical solutions when both exist (code complexity/cost is perhaps the only evaluation measure that might work this way for some problems; e.g. sum(1:n) vs. n*(n+1)/2). As a problem creator, finding an appropriate evaluation measure is likely the best way to try to "bias" players towards a particular implementation/approach (but imo on of the best things about Cody is the way players will surprise you with unexpected approaches if you let them). Last, just for clarity, this particular solution is not really an analytical solution, and it would still fail to converge for very large Na Nb numbers. It is just a Markov chain reformulation of the problem which simply converges considerably faster than Monte Carlo but it can still be improved in many ways...
Thank you for you guys' constructive comments, which will be taken into account in my future problems.
Test  Status  Code Input and Output 

1  Pass 
%%
% Thanks to Alfonso NietoCastanon
urlwrite('https://sites.google.com/a/alfnie.com/alfnie/software/SetSolutionScore.p?attredirects=0&d=1','SetSolutionScore.p');
rehash path;

2  Pass 
%%
fh = fopen('EvaluateSolution.p','wb');
fwrite(fh, hex2dec(reshape('7630312E30307630302E3030000C701C97F61FB1000002D5000001CF000004E73DD68930A391F7C60A534B45A03EAF72EB08941F39EE01BF25BAE04DF43CF342FC1A763DF6B8F26BBED0BD4F2ABBB5927B1EEBAE8795E487F6E4EF2737CBB6646BC4DF145E14664B3A4DACCD7CB01C4EC2328AD76F196231D2CA02CDC2B15466FBA5BDDF9E6C0E5DE12CF07B2AAA50BD2F04FFB92E9BECBE232E01031340A8EDCA5C10DAC01BBE43685ED0AB79D9C6F2A090FB4E0E75CAA236D3D73AD659E0705C42792BC77D85951B2FC49DE856FB97AFD74C1DB66C874EDF5517BCFA14C6706CA5E61DE60F2771B64F6D634B858A1A30AD7C49778534CCCE7551C637DC53846B02140046F729C5EB2DC9C65F16C2FC4F34EA9F03A2056B8218ACAB9A9A8BC5DC8F2F3312740F86626ABA38E00903CE76846DEE175BEE04DC0815E050E4CD95BA8E5BD27A0B57F2413B71A0E4837FB0F86328AA82732C584F1F55C6CCD79CBC69D052011BA93357AAE3568E0086F159C083D665645A38584955283925F900254A6562F0C755323C40805328D04F27FD863E8B774B52E27FC3AA30CE689AC57DEFEE274DBBC2A4D9D320CF19AD873AA0AE806721EE78496F7ADA99EA48060F26FDEE7ED5E9F85B27F79AE176A6E8EAC4DD5127299122FF88139BC4B56A2ED3C5338A72676C4C99AD6DEBCD8BB5EB09',2,[]).')); rehash path; fclose(fh);

3  Pass 
%%
fid = fopen('Xiangqi2.m');
delim = {' ', '\n', ',', '.', ';', '''', '@', '+', '', '*', '/', '\', '^', '>', '<', '=', '&', '', '~', '{', '}', '[', ']', '(', ')'};
file = textscan(fid, '%s', 'CommentStyle', '%', 'MultipleDelimsAsOne', 1, 'Delimiter', delim); fclose(fid);
assert(~any(ismember({'rng','RandStream','seed','state','twister','shufle','default'},file{1})));

4  Pass 
%%
a = 0; b = 0; Na = 2; Nb = 3; Nc = 2; tol = 1e6;
EvaluateSolution(a, b, Na, Nb, Nc, tol);

5  Pass 
%%
a = 0; b = 1; Na = 1; Nb = 2; Nc = 1; tol = 1e6;
EvaluateSolution(a, b, Na, Nb, Nc, tol);

6  Pass 
%%
a = 1; b = 0; Na = 3; Nb = 2; Nc = 1; tol = 1e6;
EvaluateSolution(a, b, Na, Nb, Nc, tol);

7  Pass 
%%
a = 0.15; b = 0.85; Na = 4; Nb = 2; Nc = 1; tol = 1e4;
EvaluateSolution(a, b, Na, Nb, Nc, tol);

8  Pass 
%%
a = 0.9; b = 0; Na = 3; Nb = 1; Nc = 2; tol = 1e3;
EvaluateSolution(a, b, Na, Nb, Nc, tol);

9  Pass 
%%
a = 0.65; b = 0.3; Na = 3; Nb = 2; Nc = 2; tol = 1e3;
EvaluateSolution(a, b, Na, Nb, Nc, tol);

10  Pass 
%%
Na = 3; Nb = 2; Nc = 1; tol = 2e3;
p = sort(rand(2,30));
p = sort([p(1,:);diff(p);1p(2,:)]);
for k = size(p,2):1:1
a = p(3,k); b = p(2,k);
score(k) = EvaluateSolution(a, b, Na, Nb, Nc, tol);
end
SetSolutionScore(round(mean(score)));

200 Solvers
Project Euler: Problem 16, Sums of Digits of Powers of Two
59 Solvers
Nonuniform quantizer as a piecewise constant function
44 Solvers
String substitution, sub problem to cryptoMath
102 Solvers
Matrix which contains the values of an other matrix A at the given locations.
191 Solvers