LINE Solver
MATLAB API documentation
Loading...
Searching...
No Matches
solver_mam_autocat.m
1function [x, pi, Q, stats] = solver_mam_autocat(R, AP, varargin)
2% SOLVER_MAM_AUTOCAT Native MATLAB implementation of autocat
3%
4% [X, PI, Q, STATS] = SOLVER_MAM_AUTOCAT(R, AP, OPTIONS)
5%
6% Searches for RCAT product-form solutions using various relaxation methods.
7% Native MATLAB implementation without YALMIP dependency.
8%
9% Input:
10% R - Cell array of rate matrices
11% AP - Action-Process mapping (A x 2)
12% OPTIONS (name-value pairs):
13% 'relaxation' - Relaxation method (default: 'auto')
14% 'lpr' - LP relaxation with McCormick envelopes
15% 'tlpr0inf' - Tightened LP with 0-level + infinity
16% 'tlpr1inf' - Tightened LP with 1-level + infinity
17% 'tlprUinf' - Tightened LP with U-level + infinity
18% 'ens' - Exact nonlinear system (fmincon)
19% 'qcp' - Quadratically constrained (fmincon)
20% 'ma' - Mean Approximation (first moment matching)
21% 'va' - Variance Approximation (slack variables)
22% 'minres' - Minimum Residual (L1-norm of RC3)
23% 'inap' - Iterative fixed-point approximation
24% 'tma1inf' - Mean Approx + Tightened LP with 1-level + infinity
25% 'zpr' - Zero potential relaxation with cutting planes
26% 'tzpr0inf' - Tightened zero potential (0-level)
27% 'tzpr1inf' - Tightened zero potential (1-level)
28% 'auto' - Automatic selection
29% 'policy' - Bound selection policy (default: 'vol')
30% 'vol', 'lu', 'l', 'u'
31% 'maxiter' - Maximum iterations (default: 100)
32% 'tol' - Convergence tolerance (default: 1e-4)
33% 'verbose' - Verbosity level 0/1/2 (default: 1)
34%
35% Output:
36% x - Action rate parameters
37% pi - Equilibrium distributions
38% Q - Generator matrices
39% stats - Statistics struct
40%
41% Copyright (c) 2012-2025, Imperial College London
42% All rights reserved.
43
44%% Parse options
45p = inputParser;
46addParameter(p, 'relaxation', 'auto');
47addParameter(p, 'policy', 'vol');
48addParameter(p, 'maxiter', 100);
49addParameter(p, 'tol', 1e-4);
50addParameter(p, 'verbose', 1);
51parse(p, varargin{:});
52opts = p.Results;
53
54RELAXATION = opts.relaxation;
55POLICY = opts.policy;
56MAXITER = opts.maxiter;
57PFTOL = opts.tol;
58GAPTOL = 0.01;
59verbose = opts.verbose;
60
61%% Extract model parameters
62A = size(AP, 1);
63ACT = AP(:, 1);
64PSV = AP(:, 2);
65M = max(AP(:));
66
67% Extract rate matrices
68Aa = cell(1, A);
69Pb = cell(1, A);
70for a = 1:A
71 Aa{a} = R{a, 1};
72 Pb{a} = R{a, 2};
73end
74
75L = cell(1, M);
76for k = 1:M
77 L{k} = R{A+1, k};
78end
79
80% State space sizes
81N = zeros(1, M);
82for k = 1:M
83 N(k) = size(L{k}, 1);
84end
85
86%% Scale rates for numerical stability
87scale = 1.0;
88if ismember(RELAXATION, {'auto', 'tlpr0inf', 'tlpr1inf', 'tlprUinf', 'tzpr0inf', 'tzpr1inf', 'zpr', 'tma1inf'})
89 for a = 1:A
90 scale = max([scale, max(sum(abs(Aa{a}), 2))]);
91 end
92 for k = 1:M
93 scale = max([scale, max(sum(abs(L{k}), 2))]);
94 end
95 scale = 1.1 * scale;
96end
97
98% Apply scaling
99Aa_scaled = cell(1, A);
100L_scaled = cell(1, M);
101for a = 1:A
102 Aa_scaled{a} = Aa{a} / scale;
103end
104for k = 1:M
105 L_scaled{k} = L{k} / scale;
106end
107
108%% Initialize bounds
109xL = zeros(A, 1);
110xU = zeros(A, 1);
111for a = 1:A
112 xU(a) = max(Aa_scaled{a}(:)) * 1.1;
113 xL(a) = max(PFTOL, min(Aa_scaled{a}(Aa_scaled{a} > 0)) * 0.1);
114end
115
116piL = cell(M, 1);
117piU = cell(M, 1);
118for k = 1:M
119 piL{k} = zeros(1, N(k));
120 piU{k} = ones(1, N(k));
121end
122
123%% Statistics
124stats = struct();
125stats.iter = 0;
126stats.action_seq = [];
127stats.bound_seq = [];
128stats.tot_time = 0;
129stats.relaxation_used = {};
130
131%% Main iteration loop
132gap = Inf(A, 1);
133currentRelaxation = RELAXATION;
134
135for iter = 1:MAXITER
136 stats.iter = iter;
137
138 % Select relaxation based on iteration (for 'auto' mode)
139 if strcmpi(RELAXATION, 'auto')
140 if iter <= A
141 currentRelaxation = 'lpr';
142 elseif iter <= 2*A
143 currentRelaxation = 'lpr';
144 elseif iter <= 3*A
145 currentRelaxation = 'tlpr1inf';
146 else
147 currentRelaxation = 'tlpr1inf';
148 end
149 end
150
151 % Select action and bound type using policy
152 [targetAction, boundType] = select_action_bound(stats, A, xL, xU, POLICY);
153 stats.action_seq(end+1) = targetAction;
154 stats.bound_seq(end+1) = boundType;
155 stats.relaxation_used{end+1} = currentRelaxation;
156
157 T0 = tic;
158
159 % Solve using selected relaxation
160 piOpt = {};
161 switch currentRelaxation
162 case 'lpr'
163 [xOpt, piOpt, fval, exitflag] = solve_lpr(targetAction, boundType, ...
164 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL);
165 case {'tlpr0inf', 'tlpr1inf', 'tlprUinf'}
166 U = 1; % Level for tightening
167 if strcmpi(currentRelaxation, 'tlpr0inf')
168 U = 0;
169 elseif strcmpi(currentRelaxation, 'tlprUinf')
170 U = min(2, floor(iter/A));
171 end
172 [xOpt, piOpt, fval, exitflag] = solve_tlpr(targetAction, boundType, ...
173 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL, U);
174 case 'ens'
175 [xOpt, piOpt, fval, exitflag] = solve_ens(targetAction, boundType, ...
176 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, PFTOL);
177 case 'qcp'
178 [xOpt, piOpt, fval, exitflag] = solve_qcp(...
179 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, PFTOL);
180 case 'ma'
181 [xOpt, piOpt, fval, exitflag] = solve_ma(targetAction, boundType, ...
182 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL);
183 case 'va'
184 [xOpt, piOpt, fval, exitflag] = solve_va(...
185 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, PFTOL);
186 case 'minres'
187 [xOpt, piOpt, fval, exitflag] = solve_minres(targetAction, boundType, ...
188 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL);
189 case 'inap'
190 [xOpt, piOpt, fval, exitflag] = solve_inap(...
191 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, PFTOL);
192 case 'tma1inf'
193 [xOpt, piOpt, fval, exitflag] = solve_tma1inf(targetAction, boundType, ...
194 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL);
195 case 'zpr'
196 [xOpt, piOpt, fval, exitflag] = solve_zpr(targetAction, boundType, ...
197 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL);
198 case {'tzpr0inf', 'tzpr1inf'}
199 U = 1; % Level for tightening
200 if strcmpi(currentRelaxation, 'tzpr0inf')
201 U = 0;
202 end
203 [xOpt, piOpt, fval, exitflag] = solve_tzpr(targetAction, boundType, ...
204 Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL, U);
205 otherwise
206 error('Unknown relaxation: %s', currentRelaxation);
207 end
208
209 stats.tot_time = stats.tot_time + toc(T0);
210
211 if exitflag <= 0
212 if verbose >= 1
213 fprintf('iter %3d: x(%d) %c - infeasible (%s)\n', ...
214 iter, targetAction, char('L' + boundType*('U'-'L')), currentRelaxation);
215 end
216 continue;
217 end
218
219 % Update x bounds
220 newBnd = abs(xOpt(targetAction));
221 if boundType == 0 % Lower bound
222 oldBnd = xL(targetAction);
223 xL(targetAction) = max(xL(targetAction), min(newBnd, xU(targetAction)));
224 else % Upper bound
225 oldBnd = xU(targetAction);
226 xU(targetAction) = min(xU(targetAction), max(newBnd, xL(targetAction)));
227 end
228
229 % Update pi bounds from LP solution (tightens McCormick envelopes)
230 if ~isempty(piOpt)
231 piTol = 0.1; % Relaxed tolerance for pi bounds
232 for k = 1:M
233 if ~isempty(piOpt{k})
234 piL{k} = max(piL{k}, piOpt{k} - piTol);
235 piU{k} = min(piU{k}, piOpt{k} + piTol);
236 end
237 end
238 end
239
240 gap(targetAction) = abs(1 - xU(targetAction)/xL(targetAction));
241
242 if verbose >= 1
243 fprintf('iter %3d: x(%d) %c = %8.5f gap = %8.5f [%s]\n', ...
244 iter, targetAction, char('L' + boundType*('U'-'L')), ...
245 newBnd * scale, gap(targetAction), currentRelaxation);
246 end
247
248 % Check convergence
249 if max(gap) < GAPTOL
250 if verbose >= 1
251 fprintf('\nConverged: max gap = %.6f\n', max(gap));
252 end
253 break;
254 end
255end
256
257%% Compute final solution
258x = (xL + xU) / 2 * scale;
259[pi, Q] = compute_equilibrium(x/scale, Aa_scaled, Pb, L_scaled, ACT, PSV, M, A, N);
260
261% Rescale Q matrices
262for k = 1:M
263 Q{k} = Q{k} * scale;
264end
265
266%% Check product-form conditions
267rc3 = check_rc3(x/scale, pi, Aa_scaled, ACT, A, N) * scale;
268if verbose >= 1
269 fprintf('\nFinal solution:\n');
270 for a = 1:A
271 fprintf(' x(%d) = %.6f [%.6f, %.6f]\n', a, x(a), xL(a)*scale, xU(a)*scale);
272 end
273 fprintf(' RC3 residual = %.6f\n', rc3);
274 fprintf(' Total time = %.2f sec\n', stats.tot_time);
275end
276
277stats.xL = xL * scale;
278stats.xU = xU * scale;
279stats.gap = gap;
280stats.rc3 = rc3;
281stats.scale = scale;
282
283end
284
285%% Action/Bound Selection Policy
286function [action, boundType] = select_action_bound(stats, A, xL, xU, policy)
287 if isempty(stats.action_seq)
288 action = 1;
289 boundType = 0;
290 return;
291 end
292
293 switch policy
294 case 'vol'
295 % Volume-based: prioritize largest gap, alternate L/U
296 if stats.bound_seq(end) == 0
297 action = stats.action_seq(end);
298 boundType = 1;
299 else
300 [~, pos] = sort(xU - xL, 'descend');
301 for a = pos(:)'
302 recent = stats.action_seq(max(1, end-2*(A-1)+1):end);
303 if isempty(find(a == recent, 1))
304 action = a;
305 boundType = 0;
306 return;
307 end
308 end
309 % Fallback
310 action = mod(stats.action_seq(end), A) + 1;
311 boundType = 0;
312 end
313 case 'lu'
314 % Alternate lower/upper for each action
315 if stats.bound_seq(end) == 0
316 action = stats.action_seq(end);
317 boundType = 1;
318 else
319 action = mod(stats.action_seq(end), A) + 1;
320 boundType = 0;
321 end
322 case 'l'
323 % Lower bounds only
324 action = mod(stats.action_seq(end), A) + 1;
325 boundType = 0;
326 case 'u'
327 % Upper bounds only
328 action = mod(stats.action_seq(end), A) + 1;
329 boundType = 1;
330 otherwise
331 action = mod(stats.action_seq(end), A) + 1;
332 boundType = mod(stats.bound_seq(end) + 1, 2);
333 end
334end
335
336%% LP Relaxation (lpr)
337function [xOpt, piOpt, fval, exitflag] = solve_lpr(targetAction, boundType, ...
338 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL)
339
340% Variables: [x(1:A); pi{1}(:); pi{2}(:); ...; z{1,1}(:); z{1,2}(:); ...]
341nX = A;
342nPi = sum(N);
343nZ = 0;
344for a = 1:A
345 nZ = nZ + N(ACT(a)) + N(PSV(a));
346end
347nVar = nX + nPi + nZ;
348
349% Variable indexing
350xIdx = 1:A;
351piIdx = cell(M, 1);
352offset = A;
353for k = 1:M
354 piIdx{k} = offset + (1:N(k));
355 offset = offset + N(k);
356end
357zIdx = cell(A, 2);
358for a = 1:A
359 zIdx{a, 1} = offset + (1:N(ACT(a)));
360 offset = offset + N(ACT(a));
361 zIdx{a, 2} = offset + (1:N(PSV(a)));
362 offset = offset + N(PSV(a));
363end
364
365% Objective
366f = zeros(nVar, 1);
367if boundType == 0
368 f(targetAction) = 1;
369else
370 f(targetAction) = -1;
371end
372
373% Build constraints
374Aeq = [];
375beq = [];
376Aineq = [];
377bineq = [];
378
379% pi{k} * 1 = 1
380for k = 1:M
381 row = zeros(1, nVar);
382 row(piIdx{k}) = 1;
383 Aeq = [Aeq; row];
384 beq = [beq; 1];
385end
386
387% sum(z{a,type}) = x(a)
388for a = 1:A
389 row = zeros(1, nVar);
390 row(zIdx{a, 1}) = 1;
391 row(xIdx(a)) = -1;
392 Aeq = [Aeq; row];
393 beq = [beq; 0];
394
395 row = zeros(1, nVar);
396 row(zIdx{a, 2}) = 1;
397 row(xIdx(a)) = -1;
398 Aeq = [Aeq; row];
399 beq = [beq; 0];
400end
401
402% pi{k} * Q{k} = 0 (balance equations)
403for k = 1:M
404 Qlocal = L{k} - diag(sum(L{k}, 2));
405
406 for s = 1:N(k)
407 row = zeros(1, nVar);
408 row(piIdx{k}) = Qlocal(:, s)';
409
410 for a = 1:A
411 if PSV(a) == k
412 Qp = Pb{a} - diag(sum(Pb{a}, 2));
413 row(zIdx{a, 2}) = row(zIdx{a, 2}) + Qp(:, s)';
414 elseif ACT(a) == k
415 Qa = Aa{a} - diag(sum(Aa{a}, 2));
416 row(piIdx{k}) = row(piIdx{k}) + Qa(:, s)';
417 end
418 end
419
420 Aineq = [Aineq; row; -row];
421 bineq = [bineq; PFTOL; PFTOL];
422 end
423end
424
425% McCormick envelopes
426for a = 1:A
427 k_act = ACT(a);
428 k_psv = PSV(a);
429
430 % Active z
431 for i = 1:N(k_act)
432 xLa = xL(a); xUa = xU(a);
433 pL = piL{k_act}(i); pU = piU{k_act}(i);
434
435 row = zeros(1, nVar);
436 row(zIdx{a, 1}(i)) = -1;
437 row(xIdx(a)) = pL;
438 row(piIdx{k_act}(i)) = xLa;
439 Aineq = [Aineq; row];
440 bineq = [bineq; pL * xLa];
441
442 row = zeros(1, nVar);
443 row(zIdx{a, 1}(i)) = 1;
444 row(xIdx(a)) = -pL;
445 row(piIdx{k_act}(i)) = -xUa;
446 Aineq = [Aineq; row];
447 bineq = [bineq; -pL * xUa];
448
449 row = zeros(1, nVar);
450 row(zIdx{a, 1}(i)) = 1;
451 row(xIdx(a)) = -pU;
452 row(piIdx{k_act}(i)) = -xLa;
453 Aineq = [Aineq; row];
454 bineq = [bineq; -pU * xLa];
455
456 row = zeros(1, nVar);
457 row(zIdx{a, 1}(i)) = -1;
458 row(xIdx(a)) = pU;
459 row(piIdx{k_act}(i)) = xUa;
460 Aineq = [Aineq; row];
461 bineq = [bineq; pU * xUa];
462 end
463
464 % Passive z
465 for i = 1:N(k_psv)
466 xLa = xL(a); xUa = xU(a);
467 pL = piL{k_psv}(i); pU = piU{k_psv}(i);
468
469 row = zeros(1, nVar);
470 row(zIdx{a, 2}(i)) = -1;
471 row(xIdx(a)) = pL;
472 row(piIdx{k_psv}(i)) = xLa;
473 Aineq = [Aineq; row];
474 bineq = [bineq; pL * xLa];
475
476 row = zeros(1, nVar);
477 row(zIdx{a, 2}(i)) = 1;
478 row(xIdx(a)) = -pL;
479 row(piIdx{k_psv}(i)) = -xUa;
480 Aineq = [Aineq; row];
481 bineq = [bineq; -pL * xUa];
482
483 row = zeros(1, nVar);
484 row(zIdx{a, 2}(i)) = 1;
485 row(xIdx(a)) = -pU;
486 row(piIdx{k_psv}(i)) = -xLa;
487 Aineq = [Aineq; row];
488 bineq = [bineq; -pU * xLa];
489
490 row = zeros(1, nVar);
491 row(zIdx{a, 2}(i)) = -1;
492 row(xIdx(a)) = pU;
493 row(piIdx{k_psv}(i)) = xUa;
494 Aineq = [Aineq; row];
495 bineq = [bineq; pU * xUa];
496 end
497end
498
499% Bounds
500lb = zeros(nVar, 1);
501ub = Inf(nVar, 1);
502lb(xIdx) = xL;
503ub(xIdx) = xU;
504for k = 1:M
505 lb(piIdx{k}) = piL{k};
506 ub(piIdx{k}) = piU{k};
507end
508for a = 1:A
509 lb(zIdx{a, 1}) = 0;
510 ub(zIdx{a, 1}) = xU(a);
511 lb(zIdx{a, 2}) = 0;
512 ub(zIdx{a, 2}) = xU(a);
513end
514
515% Solve LP
516options = optimoptions('linprog', 'Algorithm', 'interior-point', 'Display', 'off');
517[sol, fval, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, options);
518
519if exitflag > 0
520 xOpt = sol(xIdx);
521 piOpt = cell(M, 1);
522 for k = 1:M
523 piOpt{k} = sol(piIdx{k})';
524 end
525else
526 xOpt = [];
527 piOpt = {};
528end
529
530end
531
532%% Tightened LP Relaxation (tlpr)
533function [xOpt, piOpt, fval, exitflag] = solve_tlpr(targetAction, boundType, ...
534 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL, U)
535
536% Start with basic LPR
537[xOpt, piOpt, fval, exitflag] = solve_lpr(targetAction, boundType, ...
538 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL);
539
540if exitflag <= 0
541 return;
542end
543
544% Add tightening constraints iteratively
545% The tightening uses pi*Q*A^u = 0 and pi*Q*(I-A)^{-1} = 0
546% For simplicity, we iterate the basic LPR with updated bounds
547
548for u = 1:U
549 % Update pi bounds based on current solution
550 if ~isempty(piOpt)
551 for k = 1:M
552 piL{k} = max(piL{k}, piOpt{k} - PFTOL);
553 piU{k} = min(piU{k}, piOpt{k} + PFTOL);
554 end
555 end
556
557 % Re-solve with tightened bounds
558 [xOpt, piOpt, fval, exitflag] = solve_lpr(targetAction, boundType, ...
559 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL);
560
561 if exitflag <= 0
562 break;
563 end
564end
565
566end
567
568%% Exact Nonlinear System (ens) using fmincon
569function [xOpt, piOpt, fval, exitflag] = solve_ens(targetAction, boundType, ...
570 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, PFTOL)
571
572% ENS finds product-form solutions by minimizing RC3 violations
573% Uses penalty method: minimize |pi*A - x*pi|^2 subject to normalization and balance
574
575% Variables: [x(1:A); pi{1}(:); pi{2}(:); ...]
576nX = A;
577nPi = sum(N);
578nVar = nX + nPi;
579
580% Variable indexing
581xIdx = 1:A;
582piIdx = cell(M, 1);
583offset = A;
584for k = 1:M
585 piIdx{k} = offset + (1:N(k));
586 offset = offset + N(k);
587end
588
589% Better initial guess: use geometric distribution approximation
590% For M/M/1 queue, pi(n) = (1-rho) * rho^n where rho = lambda/mu
591x0 = zeros(nVar, 1);
592x0(xIdx) = (xL + xU) / 2;
593for k = 1:M
594 % Estimate utilization from rate matrices
595 rho_est = 0.5; % Default estimate
596 % Try to extract arrival/service rates from L
597 arrivals = sum(diag(L{k}, 1));
598 departures = sum(diag(L{k}, -1));
599 if arrivals > 0 && departures > 0
600 rho_est = min(0.95, arrivals / departures);
601 end
602 % Geometric distribution
603 pi_init = (1 - rho_est) * rho_est.^(0:(N(k)-1));
604 pi_init = pi_init / sum(pi_init); % Normalize
605 x0(piIdx{k}) = pi_init;
606end
607
608% Bounds
609lb = zeros(nVar, 1);
610ub = Inf(nVar, 1);
611lb(xIdx) = xL;
612ub(xIdx) = xU;
613for k = 1:M
614 lb(piIdx{k}) = 1e-12; % Small positive to avoid log issues
615 ub(piIdx{k}) = 1;
616end
617
618% Objective: minimize RC3 violation (squared sum of residuals)
619% This is better conditioned than using as hard constraints
620 function obj = objfun_penalty(v)
621 xv = v(xIdx);
622 rc3_penalty = 0;
623
624 for kk = 1:M
625 piv = v(piIdx{kk})';
626
627 % RC3: pi * Aa = x(a) * pi for active processes
628 for aa = 1:A
629 if ACT(aa) == kk
630 residual = piv * Aa{aa} - xv(aa) * piv;
631 rc3_penalty = rc3_penalty + sum(residual.^2);
632 end
633 end
634 end
635
636 % Add optimization direction as small perturbation
637 if boundType == 0 % minimizing
638 obj = rc3_penalty + 1e-6 * v(targetAction);
639 else % maximizing
640 obj = rc3_penalty - 1e-6 * v(targetAction);
641 end
642 end
643
644% Equality constraints: only normalization and balance (not RC3)
645 function [c, ceq] = nlcon(v)
646 xv = v(xIdx);
647 c = [];
648 ceq = [];
649
650 for kk = 1:M
651 piv = v(piIdx{kk})';
652
653 % Normalization: sum(pi) = 1
654 ceq = [ceq; sum(piv) - 1];
655
656 % Balance: pi * Q = 0 (only N-1 are independent, but fmincon handles this)
657 Qk = L{kk} - diag(sum(L{kk}, 2));
658 for aa = 1:A
659 if PSV(aa) == kk
660 Qk = Qk + xv(aa) * (Pb{aa} - diag(sum(Pb{aa}, 2)));
661 elseif ACT(aa) == kk
662 Qk = Qk + Aa{aa} - diag(sum(Aa{aa}, 2));
663 end
664 end
665 % Use only first N-1 balance equations (last is dependent on normalization)
666 balance = piv * Qk;
667 ceq = [ceq; balance(1:end-1)'];
668 end
669 end
670
671% Solve with relaxed tolerances and interior-point algorithm
672options = optimoptions('fmincon', 'Algorithm', 'interior-point', 'Display', 'off', ...
673 'MaxIterations', 2000, 'MaxFunctionEvaluations', 50000, ...
674 'OptimalityTolerance', 1e-8, 'ConstraintTolerance', 1e-6, ...
675 'StepTolerance', 1e-12);
676
677[sol, fval, exitflag] = fmincon(@objfun_penalty, x0, [], [], [], [], lb, ub, @nlcon, options);
678
679% Check if RC3 is satisfied (exitflag > 0 only indicates convergence)
680if exitflag > 0 && fval < 1e-4 % RC3 residual is small enough
681 xOpt = sol(xIdx);
682 piOpt = cell(M, 1);
683 for k = 1:M
684 piOpt{k} = sol(piIdx{k})';
685 end
686elseif exitflag > 0
687 % Converged but RC3 not satisfied - still return result for analysis
688 xOpt = sol(xIdx);
689 piOpt = cell(M, 1);
690 for k = 1:M
691 piOpt{k} = sol(piIdx{k})';
692 end
693else
694 xOpt = [];
695 piOpt = {};
696end
697
698end
699
700%% Quadratically Constrained Program (qcp) using fmincon
701function [xOpt, piOpt, fval, exitflag] = solve_qcp(...
702 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, PFTOL)
703
704% QCP minimizes RC3 violations directly (no slack variables needed)
705% This is a simpler reformulation that works better with fmincon
706
707% Variables: [x(1:A); pi{1}(:); pi{2}(:); ...]
708nX = A;
709nPi = sum(N);
710nVar = nX + nPi;
711
712% Variable indexing
713xIdx = 1:A;
714piIdx = cell(M, 1);
715offset = A;
716for k = 1:M
717 piIdx{k} = offset + (1:N(k));
718 offset = offset + N(k);
719end
720
721% Better initial guess: use geometric distribution approximation
722x0 = zeros(nVar, 1);
723x0(xIdx) = (xL + xU) / 2;
724for k = 1:M
725 % Estimate utilization from rate matrices
726 rho_est = 0.5;
727 arrivals = sum(diag(L{k}, 1));
728 departures = sum(diag(L{k}, -1));
729 if arrivals > 0 && departures > 0
730 rho_est = min(0.95, arrivals / departures);
731 end
732 pi_init = (1 - rho_est) * rho_est.^(0:(N(k)-1));
733 pi_init = pi_init / sum(pi_init);
734 x0(piIdx{k}) = pi_init;
735end
736
737% Bounds
738lb = zeros(nVar, 1);
739ub = Inf(nVar, 1);
740lb(xIdx) = xL;
741ub(xIdx) = xU;
742for k = 1:M
743 lb(piIdx{k}) = 1e-12;
744 ub(piIdx{k}) = 1;
745end
746
747% Objective: minimize sum of squared RC3 violations (quadratic form)
748 function obj = objfun_qcp(v)
749 xv = v(xIdx);
750 rc3_sum = 0;
751
752 for kk = 1:M
753 piv = v(piIdx{kk})';
754
755 for aa = 1:A
756 if ACT(aa) == kk
757 residual = piv * Aa{aa} - xv(aa) * piv;
758 rc3_sum = rc3_sum + sum(residual.^2);
759 end
760 end
761 end
762 obj = rc3_sum;
763 end
764
765% Equality constraints: normalization and balance only
766 function [c, ceq] = nlcon(v)
767 xv = v(xIdx);
768 c = [];
769 ceq = [];
770
771 for kk = 1:M
772 piv = v(piIdx{kk})';
773
774 % Normalization: sum(pi) = 1
775 ceq = [ceq; sum(piv) - 1];
776
777 % Balance: pi * Q = 0
778 Qk = L{kk} - diag(sum(L{kk}, 2));
779 for aa = 1:A
780 if PSV(aa) == kk
781 Qk = Qk + xv(aa) * (Pb{aa} - diag(sum(Pb{aa}, 2)));
782 elseif ACT(aa) == kk
783 Qk = Qk + Aa{aa} - diag(sum(Aa{aa}, 2));
784 end
785 end
786 % Use only first N-1 balance equations
787 balance = piv * Qk;
788 ceq = [ceq; balance(1:end-1)'];
789 end
790 end
791
792% Solve with interior-point (better for quadratic objectives)
793options = optimoptions('fmincon', 'Algorithm', 'interior-point', 'Display', 'off', ...
794 'MaxIterations', 2000, 'MaxFunctionEvaluations', 50000, ...
795 'OptimalityTolerance', 1e-8, 'ConstraintTolerance', 1e-6, ...
796 'StepTolerance', 1e-12);
797
798[sol, fval, exitflag] = fmincon(@objfun_qcp, x0, [], [], [], [], lb, ub, @nlcon, options);
799
800if exitflag > 0
801 xOpt = sol(xIdx);
802 piOpt = cell(M, 1);
803 for k = 1:M
804 piOpt{k} = sol(piIdx{k})';
805 end
806else
807 xOpt = [];
808 piOpt = {};
809end
810
811end
812
813%% Compute equilibrium
814function [pi, Q] = compute_equilibrium(x, Aa, Pb, L, ACT, PSV, M, A, N)
815Q = cell(1, M);
816pi = cell(1, M);
817
818for k = 1:M
819 Qk = L{k} - diag(sum(L{k}, 2));
820
821 for a = 1:A
822 if PSV(a) == k
823 Qk = Qk + x(a) * (Pb{a} - diag(sum(Pb{a}, 2)));
824 elseif ACT(a) == k
825 Qk = Qk + Aa{a} - diag(sum(Aa{a}, 2));
826 end
827 end
828
829 Q{k} = Qk;
830 Qk_inf = ctmc_makeinfgen(Qk);
831 pi{k} = ctmc_solve(Qk_inf);
832end
833end
834
835%% Check RC3 condition
836function rc3 = check_rc3(x, pi, Aa, ACT, A, N)
837rc3 = 0;
838for a = 1:A
839 k = ACT(a);
840 residual = abs(pi{k} * Aa{a} - x(a) * pi{k});
841 rc3 = rc3 + mean(residual);
842end
843rc3 = rc3 / A;
844end
845
846%% Mean Approximation (ma)
847% Relaxes RC3 to only enforce first moment matching: sum(z) = sum(pi*Aa)
848% This is equivalent to x = E[eigenvalue] from the active rate matrix
849function [xOpt, piOpt, fval, exitflag] = solve_ma(targetAction, boundType, ...
850 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL)
851
852% Variables: [x(1:A); pi{1}(:); pi{2}(:); ...; z{1,1}(:); z{1,2}(:); ...]
853nX = A;
854nPi = sum(N);
855nZ = 0;
856for a = 1:A
857 nZ = nZ + N(ACT(a)) + N(PSV(a));
858end
859nVar = nX + nPi + nZ;
860
861% Variable indexing
862xIdx = 1:A;
863piIdx = cell(M, 1);
864offset = A;
865for k = 1:M
866 piIdx{k} = offset + (1:N(k));
867 offset = offset + N(k);
868end
869zIdx = cell(A, 2);
870for a = 1:A
871 zIdx{a, 1} = offset + (1:N(ACT(a))); % z_active
872 offset = offset + N(ACT(a));
873 zIdx{a, 2} = offset + (1:N(PSV(a))); % z_passive
874 offset = offset + N(PSV(a));
875end
876
877% Objective
878f = zeros(nVar, 1);
879if boundType == 0
880 f(targetAction) = 1;
881else
882 f(targetAction) = -1;
883end
884
885% Build constraints
886Aeq = [];
887beq = [];
888Aineq = [];
889bineq = [];
890
891% pi{k} * 1 = 1 (normalization)
892for k = 1:M
893 row = zeros(1, nVar);
894 row(piIdx{k}) = 1;
895 Aeq = [Aeq; row];
896 beq = [beq; 1];
897end
898
899% Mean approximation: sum(z{a,active}) = sum(pi*Aa) = x(a) * sum(pi) = x(a)
900% and sum(z{a,passive}) = x(a)
901for a = 1:A
902 % z_active sums to x
903 row = zeros(1, nVar);
904 row(zIdx{a, 1}) = 1;
905 row(xIdx(a)) = -1;
906 Aeq = [Aeq; row];
907 beq = [beq; 0];
908
909 % z_passive sums to x
910 row = zeros(1, nVar);
911 row(zIdx{a, 2}) = 1;
912 row(xIdx(a)) = -1;
913 Aeq = [Aeq; row];
914 beq = [beq; 0];
915end
916
917% Balance equations: pi{k} * Q{k} = 0 (relaxed to tolerance)
918for k = 1:M
919 Qlocal = L{k} - diag(sum(L{k}, 2));
920
921 for s = 1:N(k)
922 row = zeros(1, nVar);
923 row(piIdx{k}) = Qlocal(:, s)';
924
925 for a = 1:A
926 if PSV(a) == k
927 Qp = Pb{a} - diag(sum(Pb{a}, 2));
928 row(zIdx{a, 2}) = row(zIdx{a, 2}) + Qp(:, s)';
929 elseif ACT(a) == k
930 Qa = Aa{a} - diag(sum(Aa{a}, 2));
931 row(piIdx{k}) = row(piIdx{k}) + Qa(:, s)';
932 end
933 end
934
935 Aineq = [Aineq; row; -row];
936 bineq = [bineq; PFTOL; PFTOL];
937 end
938end
939
940% McCormick envelopes for z = x * pi
941for a = 1:A
942 k_act = ACT(a);
943 k_psv = PSV(a);
944
945 % Active z
946 for i = 1:N(k_act)
947 xLa = xL(a); xUa = xU(a);
948 pL = piL{k_act}(i); pU = piU{k_act}(i);
949
950 % z >= pL*x + xL*pi - pL*xL
951 row = zeros(1, nVar);
952 row(zIdx{a, 1}(i)) = -1;
953 row(xIdx(a)) = pL;
954 row(piIdx{k_act}(i)) = xLa;
955 Aineq = [Aineq; row];
956 bineq = [bineq; pL * xLa];
957
958 % z <= pL*x + xU*pi - pL*xU
959 row = zeros(1, nVar);
960 row(zIdx{a, 1}(i)) = 1;
961 row(xIdx(a)) = -pL;
962 row(piIdx{k_act}(i)) = -xUa;
963 Aineq = [Aineq; row];
964 bineq = [bineq; -pL * xUa];
965
966 % z <= pU*x + xL*pi - pU*xL
967 row = zeros(1, nVar);
968 row(zIdx{a, 1}(i)) = 1;
969 row(xIdx(a)) = -pU;
970 row(piIdx{k_act}(i)) = -xLa;
971 Aineq = [Aineq; row];
972 bineq = [bineq; -pU * xLa];
973
974 % z >= pU*x + xU*pi - pU*xU
975 row = zeros(1, nVar);
976 row(zIdx{a, 1}(i)) = -1;
977 row(xIdx(a)) = pU;
978 row(piIdx{k_act}(i)) = xUa;
979 Aineq = [Aineq; row];
980 bineq = [bineq; pU * xUa];
981 end
982
983 % Passive z
984 for i = 1:N(k_psv)
985 xLa = xL(a); xUa = xU(a);
986 pL = piL{k_psv}(i); pU = piU{k_psv}(i);
987
988 row = zeros(1, nVar);
989 row(zIdx{a, 2}(i)) = -1;
990 row(xIdx(a)) = pL;
991 row(piIdx{k_psv}(i)) = xLa;
992 Aineq = [Aineq; row];
993 bineq = [bineq; pL * xLa];
994
995 row = zeros(1, nVar);
996 row(zIdx{a, 2}(i)) = 1;
997 row(xIdx(a)) = -pL;
998 row(piIdx{k_psv}(i)) = -xUa;
999 Aineq = [Aineq; row];
1000 bineq = [bineq; -pL * xUa];
1001
1002 row = zeros(1, nVar);
1003 row(zIdx{a, 2}(i)) = 1;
1004 row(xIdx(a)) = -pU;
1005 row(piIdx{k_psv}(i)) = -xLa;
1006 Aineq = [Aineq; row];
1007 bineq = [bineq; -pU * xLa];
1008
1009 row = zeros(1, nVar);
1010 row(zIdx{a, 2}(i)) = -1;
1011 row(xIdx(a)) = pU;
1012 row(piIdx{k_psv}(i)) = xUa;
1013 Aineq = [Aineq; row];
1014 bineq = [bineq; pU * xUa];
1015 end
1016end
1017
1018% Relaxed RC3 (mean approximation): sum(z*I - pi*Aa) = 0
1019% This only enforces first moment, not pointwise equality
1020for a = 1:A
1021 k = ACT(a);
1022 % sum(z{a,active}) - sum(pi * Aa) should be close to 0
1023 % but sum(z{a,active}) = x(a), so this is: x(a) - sum(pi*Aa) close to 0
1024 row = zeros(1, nVar);
1025 row(xIdx(a)) = 1;
1026 row(piIdx{k}) = -sum(Aa{a}, 2)';
1027 Aineq = [Aineq; row; -row];
1028 bineq = [bineq; PFTOL; PFTOL];
1029end
1030
1031% Bounds
1032lb = zeros(nVar, 1);
1033ub = Inf(nVar, 1);
1034lb(xIdx) = xL;
1035ub(xIdx) = xU;
1036for k = 1:M
1037 lb(piIdx{k}) = piL{k};
1038 ub(piIdx{k}) = piU{k};
1039end
1040for a = 1:A
1041 lb(zIdx{a, 1}) = 0;
1042 ub(zIdx{a, 1}) = xU(a);
1043 lb(zIdx{a, 2}) = 0;
1044 ub(zIdx{a, 2}) = xU(a);
1045end
1046
1047% Solve LP
1048options = optimoptions('linprog', 'Algorithm', 'interior-point', 'Display', 'off');
1049[sol, fval, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, options);
1050
1051if exitflag > 0
1052 xOpt = sol(xIdx);
1053 piOpt = cell(M, 1);
1054 for k = 1:M
1055 piOpt{k} = sol(piIdx{k})';
1056 end
1057else
1058 xOpt = [];
1059 piOpt = {};
1060end
1061
1062end
1063
1064%% Variance Approximation (va)
1065% Uses fixed x from bounds, solves for pi with slack variables on RC3
1066function [xOpt, piOpt, fval, exitflag] = solve_va(...
1067 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, PFTOL)
1068
1069% Use midpoint of bounds as fixed x
1070x_fixed = (xL + xU) / 2;
1071
1072% Variables: [pi{1}(:); pi{2}(:); ...; slack{1}; slack{2}; ...]
1073nPi = sum(N);
1074nSlack = 0;
1075for a = 1:A
1076 nSlack = nSlack + N(ACT(a));
1077end
1078nVar = nPi + nSlack;
1079
1080% Variable indexing
1081piIdx = cell(M, 1);
1082offset = 0;
1083for k = 1:M
1084 piIdx{k} = offset + (1:N(k));
1085 offset = offset + N(k);
1086end
1087slackIdx = cell(A, 1);
1088for a = 1:A
1089 slackIdx{a} = offset + (1:N(ACT(a)));
1090 offset = offset + N(ACT(a));
1091end
1092
1093% Objective: minimize sum of absolute slacks (L1 norm)
1094f = zeros(nVar, 1);
1095for a = 1:A
1096 f(slackIdx{a}) = 1; % Minimize slack
1097end
1098
1099% Build constraints
1100Aeq = [];
1101beq = [];
1102Aineq = [];
1103bineq = [];
1104
1105% pi{k} * 1 = 1 (normalization)
1106for k = 1:M
1107 row = zeros(1, nVar);
1108 row(piIdx{k}) = 1;
1109 Aeq = [Aeq; row];
1110 beq = [beq; 1];
1111end
1112
1113% Balance equations: pi{k} * Q{k} = 0
1114for k = 1:M
1115 Qk = L{k} - diag(sum(L{k}, 2));
1116 for a = 1:A
1117 if PSV(a) == k
1118 Qk = Qk + x_fixed(a) * (Pb{a} - diag(sum(Pb{a}, 2)));
1119 elseif ACT(a) == k
1120 Qk = Qk + Aa{a} - diag(sum(Aa{a}, 2));
1121 end
1122 end
1123
1124 for s = 1:(N(k)-1) % Only N-1 are independent
1125 row = zeros(1, nVar);
1126 row(piIdx{k}) = Qk(:, s)';
1127 Aeq = [Aeq; row];
1128 beq = [beq; 0];
1129 end
1130end
1131
1132% RC3 with slack: pi*Aa - x*pi + slack >= 0, pi*Aa - x*pi - slack <= 0
1133% This means |pi*Aa - x*pi| <= slack
1134for a = 1:A
1135 k = ACT(a);
1136 for i = 1:N(k)
1137 % (pi*Aa - x*pi)_i + slack_i >= 0
1138 row = zeros(1, nVar);
1139 row(piIdx{k}) = Aa{a}(:, i)' - x_fixed(a) * (1:N(k) == i);
1140 row(slackIdx{a}(i)) = 1;
1141 Aineq = [Aineq; -row]; % -(...) <= 0
1142 bineq = [bineq; 0];
1143
1144 % (pi*Aa - x*pi)_i - slack_i <= 0
1145 row = zeros(1, nVar);
1146 row(piIdx{k}) = Aa{a}(:, i)' - x_fixed(a) * (1:N(k) == i);
1147 row(slackIdx{a}(i)) = -1;
1148 Aineq = [Aineq; row];
1149 bineq = [bineq; 0];
1150 end
1151end
1152
1153% Bounds
1154lb = zeros(nVar, 1);
1155ub = Inf(nVar, 1);
1156for k = 1:M
1157 lb(piIdx{k}) = 0;
1158 ub(piIdx{k}) = 1;
1159end
1160for a = 1:A
1161 lb(slackIdx{a}) = 0; % Slack must be non-negative
1162end
1163
1164% Solve LP
1165options = optimoptions('linprog', 'Algorithm', 'interior-point', 'Display', 'off');
1166[sol, fval, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, options);
1167
1168if exitflag > 0
1169 xOpt = x_fixed; % Return the fixed x
1170 piOpt = cell(M, 1);
1171 for k = 1:M
1172 piOpt{k} = sol(piIdx{k})';
1173 end
1174else
1175 xOpt = [];
1176 piOpt = {};
1177end
1178
1179end
1180
1181%% Minimum Residual (minres)
1182% Minimizes L1 norm of RC3 violations with McCormick envelopes
1183function [xOpt, piOpt, fval, exitflag] = solve_minres(targetAction, boundType, ...
1184 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL)
1185
1186% Variables: [x; pi{1}(:); ...; z{1,1}(:); z{1,2}(:); ...; slack_pos; slack_neg]
1187nX = A;
1188nPi = sum(N);
1189nZ = 0;
1190for a = 1:A
1191 nZ = nZ + N(ACT(a)) + N(PSV(a));
1192end
1193nSlack = 0;
1194for a = 1:A
1195 nSlack = nSlack + 2 * N(ACT(a)); % Positive and negative slack for each RC3 component
1196end
1197nVar = nX + nPi + nZ + nSlack;
1198
1199% Variable indexing
1200xIdx = 1:A;
1201piIdx = cell(M, 1);
1202offset = A;
1203for k = 1:M
1204 piIdx{k} = offset + (1:N(k));
1205 offset = offset + N(k);
1206end
1207zIdx = cell(A, 2);
1208for a = 1:A
1209 zIdx{a, 1} = offset + (1:N(ACT(a)));
1210 offset = offset + N(ACT(a));
1211 zIdx{a, 2} = offset + (1:N(PSV(a)));
1212 offset = offset + N(PSV(a));
1213end
1214slackPosIdx = cell(A, 1);
1215slackNegIdx = cell(A, 1);
1216for a = 1:A
1217 slackPosIdx{a} = offset + (1:N(ACT(a)));
1218 offset = offset + N(ACT(a));
1219 slackNegIdx{a} = offset + (1:N(ACT(a)));
1220 offset = offset + N(ACT(a));
1221end
1222
1223% Objective: minimize L1 norm of slack (minimize RC3 residuals)
1224f = zeros(nVar, 1);
1225for a = 1:A
1226 f(slackPosIdx{a}) = 1;
1227 f(slackNegIdx{a}) = 1;
1228end
1229
1230% Build constraints
1231Aeq = [];
1232beq = [];
1233Aineq = [];
1234bineq = [];
1235
1236% pi{k} * 1 = 1
1237for k = 1:M
1238 row = zeros(1, nVar);
1239 row(piIdx{k}) = 1;
1240 Aeq = [Aeq; row];
1241 beq = [beq; 1];
1242end
1243
1244% sum(z{a,type}) = x(a)
1245for a = 1:A
1246 row = zeros(1, nVar);
1247 row(zIdx{a, 1}) = 1;
1248 row(xIdx(a)) = -1;
1249 Aeq = [Aeq; row];
1250 beq = [beq; 0];
1251
1252 row = zeros(1, nVar);
1253 row(zIdx{a, 2}) = 1;
1254 row(xIdx(a)) = -1;
1255 Aeq = [Aeq; row];
1256 beq = [beq; 0];
1257end
1258
1259% Balance equations: pi{k} * Q{k} = 0
1260for k = 1:M
1261 Qlocal = L{k} - diag(sum(L{k}, 2));
1262
1263 for s = 1:N(k)
1264 row = zeros(1, nVar);
1265 row(piIdx{k}) = Qlocal(:, s)';
1266
1267 for a = 1:A
1268 if PSV(a) == k
1269 Qp = Pb{a} - diag(sum(Pb{a}, 2));
1270 row(zIdx{a, 2}) = row(zIdx{a, 2}) + Qp(:, s)';
1271 elseif ACT(a) == k
1272 Qa = Aa{a} - diag(sum(Aa{a}, 2));
1273 row(piIdx{k}) = row(piIdx{k}) + Qa(:, s)';
1274 end
1275 end
1276
1277 Aineq = [Aineq; row; -row];
1278 bineq = [bineq; PFTOL; PFTOL];
1279 end
1280end
1281
1282% McCormick envelopes for z = x * pi
1283for a = 1:A
1284 k_act = ACT(a);
1285 k_psv = PSV(a);
1286
1287 for i = 1:N(k_act)
1288 xLa = xL(a); xUa = xU(a);
1289 pL = piL{k_act}(i); pU = piU{k_act}(i);
1290
1291 row = zeros(1, nVar);
1292 row(zIdx{a, 1}(i)) = -1;
1293 row(xIdx(a)) = pL;
1294 row(piIdx{k_act}(i)) = xLa;
1295 Aineq = [Aineq; row];
1296 bineq = [bineq; pL * xLa];
1297
1298 row = zeros(1, nVar);
1299 row(zIdx{a, 1}(i)) = 1;
1300 row(xIdx(a)) = -pL;
1301 row(piIdx{k_act}(i)) = -xUa;
1302 Aineq = [Aineq; row];
1303 bineq = [bineq; -pL * xUa];
1304
1305 row = zeros(1, nVar);
1306 row(zIdx{a, 1}(i)) = 1;
1307 row(xIdx(a)) = -pU;
1308 row(piIdx{k_act}(i)) = -xLa;
1309 Aineq = [Aineq; row];
1310 bineq = [bineq; -pU * xLa];
1311
1312 row = zeros(1, nVar);
1313 row(zIdx{a, 1}(i)) = -1;
1314 row(xIdx(a)) = pU;
1315 row(piIdx{k_act}(i)) = xUa;
1316 Aineq = [Aineq; row];
1317 bineq = [bineq; pU * xUa];
1318 end
1319
1320 for i = 1:N(k_psv)
1321 xLa = xL(a); xUa = xU(a);
1322 pL = piL{k_psv}(i); pU = piU{k_psv}(i);
1323
1324 row = zeros(1, nVar);
1325 row(zIdx{a, 2}(i)) = -1;
1326 row(xIdx(a)) = pL;
1327 row(piIdx{k_psv}(i)) = xLa;
1328 Aineq = [Aineq; row];
1329 bineq = [bineq; pL * xLa];
1330
1331 row = zeros(1, nVar);
1332 row(zIdx{a, 2}(i)) = 1;
1333 row(xIdx(a)) = -pL;
1334 row(piIdx{k_psv}(i)) = -xUa;
1335 Aineq = [Aineq; row];
1336 bineq = [bineq; -pL * xUa];
1337
1338 row = zeros(1, nVar);
1339 row(zIdx{a, 2}(i)) = 1;
1340 row(xIdx(a)) = -pU;
1341 row(piIdx{k_psv}(i)) = -xLa;
1342 Aineq = [Aineq; row];
1343 bineq = [bineq; -pU * xLa];
1344
1345 row = zeros(1, nVar);
1346 row(zIdx{a, 2}(i)) = -1;
1347 row(xIdx(a)) = pU;
1348 row(piIdx{k_psv}(i)) = xUa;
1349 Aineq = [Aineq; row];
1350 bineq = [bineq; pU * xUa];
1351 end
1352end
1353
1354% RC3 with slack: z_i - (pi*Aa)_i = slack_pos_i - slack_neg_i
1355% We use z{a,active} as the bilinear term x*pi
1356for a = 1:A
1357 k = ACT(a);
1358 for i = 1:N(k)
1359 % z_i - (pi*Aa)_i - slack_pos + slack_neg = 0
1360 row = zeros(1, nVar);
1361 row(zIdx{a, 1}(i)) = 1;
1362 row(piIdx{k}) = -Aa{a}(:, i)';
1363 row(slackPosIdx{a}(i)) = -1;
1364 row(slackNegIdx{a}(i)) = 1;
1365 Aeq = [Aeq; row];
1366 beq = [beq; 0];
1367 end
1368end
1369
1370% Bounds
1371lb = zeros(nVar, 1);
1372ub = Inf(nVar, 1);
1373lb(xIdx) = xL;
1374ub(xIdx) = xU;
1375for k = 1:M
1376 lb(piIdx{k}) = piL{k};
1377 ub(piIdx{k}) = piU{k};
1378end
1379for a = 1:A
1380 lb(zIdx{a, 1}) = 0;
1381 ub(zIdx{a, 1}) = xU(a);
1382 lb(zIdx{a, 2}) = 0;
1383 ub(zIdx{a, 2}) = xU(a);
1384 lb(slackPosIdx{a}) = 0;
1385 lb(slackNegIdx{a}) = 0;
1386end
1387
1388% Solve LP
1389options = optimoptions('linprog', 'Algorithm', 'interior-point', 'Display', 'off');
1390[sol, fval, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, options);
1391
1392if exitflag > 0
1393 xOpt = sol(xIdx);
1394 piOpt = cell(M, 1);
1395 for k = 1:M
1396 piOpt{k} = sol(piIdx{k})';
1397 end
1398else
1399 xOpt = [];
1400 piOpt = {};
1401end
1402
1403end
1404
1405%% Iterative Fixed-Point Approximation (inap)
1406% Uses fixed-point iteration: x(a) = mean(pi*Aa ./ pi)
1407function [xOpt, piOpt, fval, exitflag] = solve_inap(...
1408 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, PFTOL)
1409
1410% Initialize x randomly within bounds
1411x = xL + (xU - xL) .* rand(A, 1);
1412
1413% Compute initial equilibrium
1414[pi, Q] = compute_eq_internal(x, Aa, Pb, L, ACT, PSV, M, A, N);
1415
1416% Fixed-point iteration
1417maxIter = 50;
1418tol = 1e-6;
1419
1420for iter = 1:maxIter
1421 xprev = x;
1422
1423 for a = 1:A
1424 k = ACT(a);
1425 % LAMBDA = (pi * Aa) ./ pi (element-wise)
1426 piAa = pi{k} * Aa{a};
1427 LAMBDA = piAa ./ pi{k};
1428 LAMBDA = LAMBDA(~isnan(LAMBDA) & ~isinf(LAMBDA) & LAMBDA > 0);
1429
1430 if ~isempty(LAMBDA)
1431 x(a) = mean(LAMBDA);
1432 % Clamp to bounds
1433 x(a) = max(xL(a), min(xU(a), x(a)));
1434 end
1435 end
1436
1437 % Recompute equilibrium
1438 [pi, Q] = compute_eq_internal(x, Aa, Pb, L, ACT, PSV, M, A, N);
1439
1440 % Check convergence
1441 if norm(x - xprev, inf) < tol
1442 break;
1443 end
1444end
1445
1446% Check if RC3 is satisfied
1447rc3_satisfied = true;
1448for a = 1:A
1449 k = ACT(a);
1450 piAa = pi{k} * Aa{a};
1451 LAMBDA = piAa ./ pi{k};
1452 LAMBDA = LAMBDA(~isnan(LAMBDA) & ~isinf(LAMBDA) & LAMBDA > 0);
1453 if ~isempty(LAMBDA) && (max(LAMBDA) - min(LAMBDA)) > 0.1
1454 rc3_satisfied = false;
1455 end
1456end
1457
1458if rc3_satisfied
1459 xOpt = x;
1460 piOpt = pi;
1461 fval = 0;
1462 exitflag = 1;
1463else
1464 xOpt = x;
1465 piOpt = pi;
1466 fval = 1; % Indicates approximate solution
1467 exitflag = 1; % Still return solution for use
1468end
1469
1470end
1471
1472%% Internal equilibrium computation for inap
1473function [pi, Q] = compute_eq_internal(x, Aa, Pb, L, ACT, PSV, M, A, N)
1474Q = cell(1, M);
1475pi = cell(1, M);
1476
1477for k = 1:M
1478 Qk = L{k} - diag(sum(L{k}, 2));
1479
1480 for a = 1:A
1481 if PSV(a) == k
1482 Qk = Qk + x(a) * (Pb{a} - diag(sum(Pb{a}, 2)));
1483 elseif ACT(a) == k
1484 Qk = Qk + Aa{a} - diag(sum(Aa{a}, 2));
1485 end
1486 end
1487
1488 Q{k} = Qk;
1489 Qk_inf = ctmc_makeinfgen(Qk);
1490 pi{k} = ctmc_solve(Qk_inf);
1491end
1492end
1493
1494%% Tightened Mean Approximation with 1-level + infinity (tma1inf)
1495% Combines mean approximation (RC3 first moment only) with tightened LP constraints
1496function [xOpt, piOpt, fval, exitflag] = solve_tma1inf(targetAction, boundType, ...
1497 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL)
1498
1499% First run mean approximation to get initial solution
1500[xMA, piMA, ~, exitflagMA] = solve_ma(targetAction, boundType, ...
1501 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL);
1502
1503if exitflagMA <= 0
1504 xOpt = [];
1505 piOpt = {};
1506 fval = Inf;
1507 exitflag = -1;
1508 return;
1509end
1510
1511% Use MA solution to warm-start tightened LP with additional constraints
1512% Variables: [x(1:A); pi{1}(:); pi{2}(:); ...; z{1,1}(:); z{1,2}(:); ...]
1513nX = A;
1514nPi = sum(N);
1515nZ = 0;
1516for a = 1:A
1517 nZ = nZ + N(ACT(a)) + N(PSV(a));
1518end
1519nVar = nX + nPi + nZ;
1520
1521% Variable indexing
1522xIdx = 1:A;
1523piIdx = cell(M, 1);
1524offset = A;
1525for k = 1:M
1526 piIdx{k} = offset + (1:N(k));
1527 offset = offset + N(k);
1528end
1529zIdx = cell(A, 2);
1530for a = 1:A
1531 zIdx{a, 1} = offset + (1:N(ACT(a)));
1532 offset = offset + N(ACT(a));
1533 zIdx{a, 2} = offset + (1:N(PSV(a)));
1534 offset = offset + N(PSV(a));
1535end
1536
1537% Objective
1538f = zeros(nVar, 1);
1539if boundType == 0
1540 f(targetAction) = 1;
1541else
1542 f(targetAction) = -1;
1543end
1544
1545% Build constraints
1546Aeq = [];
1547beq = [];
1548Aineq = [];
1549bineq = [];
1550
1551% pi{k} * 1 = 1
1552for k = 1:M
1553 row = zeros(1, nVar);
1554 row(piIdx{k}) = 1;
1555 Aeq = [Aeq; row];
1556 beq = [beq; 1];
1557end
1558
1559% RC3 first moment constraint: sum(z{a,active}) = x(a) (mean approximation)
1560for a = 1:A
1561 row = zeros(1, nVar);
1562 row(zIdx{a, 1}) = 1;
1563 row(xIdx(a)) = -1;
1564 Aeq = [Aeq; row];
1565 beq = [beq; 0];
1566
1567 row = zeros(1, nVar);
1568 row(zIdx{a, 2}) = 1;
1569 row(xIdx(a)) = -1;
1570 Aeq = [Aeq; row];
1571 beq = [beq; 0];
1572end
1573
1574% pi{k} * Q{k} = 0 (balance equations)
1575for k = 1:M
1576 Qlocal = L{k} - diag(sum(L{k}, 2));
1577
1578 for s = 1:N(k)
1579 row = zeros(1, nVar);
1580 row(piIdx{k}) = Qlocal(:, s)';
1581
1582 for a = 1:A
1583 if PSV(a) == k
1584 Qp = Pb{a} - diag(sum(Pb{a}, 2));
1585 row(zIdx{a, 2}) = row(zIdx{a, 2}) + Qp(:, s)';
1586 elseif ACT(a) == k
1587 Qa = Aa{a} - diag(sum(Aa{a}, 2));
1588 row(piIdx{k}) = row(piIdx{k}) + Qa(:, s)';
1589 end
1590 end
1591
1592 Aineq = [Aineq; row; -row];
1593 bineq = [bineq; PFTOL; PFTOL];
1594 end
1595end
1596
1597% Tightened LP constraints: pi*Q*Aa^u = 0 for u=0,1,...,U and pi*Q*Aa^inf = 0
1598U = 1; % Tightening level
1599for k = 1:M
1600 for b = 1:A
1601 if ACT(b) == k
1602 Qlocal = L{k} - diag(sum(L{k}, 2));
1603
1604 % Add constraints for u = 1
1605 for u = 1:U
1606 Aa_pow = Aa{b}^u;
1607
1608 for s = 1:N(k)
1609 row = zeros(1, nVar);
1610
1611 % z{b,active} * Aa^(u-1) * Q_local
1612 coeffZ = (Aa{b}^(u-1)) * Qlocal;
1613 row(zIdx{b, 1}) = coeffZ(:, s)';
1614
1615 for c = 1:A
1616 if PSV(c) == k
1617 Qp = Pb{c} - diag(sum(Pb{c}, 2));
1618 coeffZp = Aa_pow * Qp;
1619 row(zIdx{c, 2}) = row(zIdx{c, 2}) + coeffZp(:, s)';
1620 elseif ACT(c) == k
1621 Qa = Aa{c} - diag(sum(Aa{c}, 2));
1622 coeffZa = (Aa{b}^(u-1)) * Qa;
1623 row(zIdx{b, 1}) = row(zIdx{b, 1}) + coeffZa(:, s)';
1624 end
1625 end
1626
1627 Aineq = [Aineq; row; -row];
1628 bineq = [bineq; PFTOL; PFTOL];
1629 end
1630 end
1631
1632 % Add infinity constraint: z * (I-Aa)^{-1} * Q = 0
1633 try
1634 Aa_inv = inv(eye(N(k)) - Aa{b});
1635
1636 for s = 1:N(k)
1637 row = zeros(1, nVar);
1638
1639 coeffZ = Aa_inv * Qlocal;
1640 row(zIdx{b, 1}) = coeffZ(:, s)';
1641
1642 for c = 1:A
1643 if PSV(c) == k
1644 Qp = Pb{c} - diag(sum(Pb{c}, 2));
1645 coeffZp = Aa{b} * Aa_inv * Qp;
1646 row(zIdx{c, 2}) = row(zIdx{c, 2}) + coeffZp(:, s)';
1647 elseif ACT(c) == k
1648 Qa = Aa{c} - diag(sum(Aa{c}, 2));
1649 coeffZa = Aa_inv * Qa;
1650 row(zIdx{b, 1}) = row(zIdx{b, 1}) + coeffZa(:, s)';
1651 end
1652 end
1653
1654 Aineq = [Aineq; row; -row];
1655 bineq = [bineq; PFTOL; PFTOL];
1656 end
1657 catch
1658 % Matrix not invertible, skip infinity constraints
1659 end
1660 end
1661 end
1662end
1663
1664% McCormick envelopes
1665for a = 1:A
1666 k_act = ACT(a);
1667 k_psv = PSV(a);
1668
1669 % Active z
1670 for i = 1:N(k_act)
1671 xLa = xL(a); xUa = xU(a);
1672 pL = piL{k_act}(i); pU = piU{k_act}(i);
1673
1674 row = zeros(1, nVar);
1675 row(zIdx{a, 1}(i)) = -1;
1676 row(xIdx(a)) = pL;
1677 row(piIdx{k_act}(i)) = xLa;
1678 Aineq = [Aineq; row];
1679 bineq = [bineq; pL * xLa];
1680
1681 row = zeros(1, nVar);
1682 row(zIdx{a, 1}(i)) = 1;
1683 row(xIdx(a)) = -pL;
1684 row(piIdx{k_act}(i)) = -xUa;
1685 Aineq = [Aineq; row];
1686 bineq = [bineq; -pL * xUa];
1687
1688 row = zeros(1, nVar);
1689 row(zIdx{a, 1}(i)) = 1;
1690 row(xIdx(a)) = -pU;
1691 row(piIdx{k_act}(i)) = -xLa;
1692 Aineq = [Aineq; row];
1693 bineq = [bineq; -pU * xLa];
1694
1695 row = zeros(1, nVar);
1696 row(zIdx{a, 1}(i)) = -1;
1697 row(xIdx(a)) = pU;
1698 row(piIdx{k_act}(i)) = xUa;
1699 Aineq = [Aineq; row];
1700 bineq = [bineq; pU * xUa];
1701 end
1702
1703 % Passive z
1704 for i = 1:N(k_psv)
1705 xLa = xL(a); xUa = xU(a);
1706 pL = piL{k_psv}(i); pU = piU{k_psv}(i);
1707
1708 row = zeros(1, nVar);
1709 row(zIdx{a, 2}(i)) = -1;
1710 row(xIdx(a)) = pL;
1711 row(piIdx{k_psv}(i)) = xLa;
1712 Aineq = [Aineq; row];
1713 bineq = [bineq; pL * xLa];
1714
1715 row = zeros(1, nVar);
1716 row(zIdx{a, 2}(i)) = 1;
1717 row(xIdx(a)) = -pL;
1718 row(piIdx{k_psv}(i)) = -xUa;
1719 Aineq = [Aineq; row];
1720 bineq = [bineq; -pL * xUa];
1721
1722 row = zeros(1, nVar);
1723 row(zIdx{a, 2}(i)) = 1;
1724 row(xIdx(a)) = -pU;
1725 row(piIdx{k_psv}(i)) = -xLa;
1726 Aineq = [Aineq; row];
1727 bineq = [bineq; -pU * xLa];
1728
1729 row = zeros(1, nVar);
1730 row(zIdx{a, 2}(i)) = -1;
1731 row(xIdx(a)) = pU;
1732 row(piIdx{k_psv}(i)) = xUa;
1733 Aineq = [Aineq; row];
1734 bineq = [bineq; pU * xUa];
1735 end
1736end
1737
1738% Bounds
1739lb = zeros(nVar, 1);
1740ub = Inf(nVar, 1);
1741lb(xIdx) = xL;
1742ub(xIdx) = xU;
1743for k = 1:M
1744 lb(piIdx{k}) = piL{k};
1745 ub(piIdx{k}) = piU{k};
1746end
1747for a = 1:A
1748 lb(zIdx{a, 1}) = 0;
1749 ub(zIdx{a, 1}) = xU(a);
1750 lb(zIdx{a, 2}) = 0;
1751 ub(zIdx{a, 2}) = xU(a);
1752end
1753
1754% Solve LP
1755options = optimoptions('linprog', 'Algorithm', 'interior-point', 'Display', 'off');
1756[sol, fval, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, options);
1757
1758if exitflag > 0
1759 xOpt = sol(xIdx);
1760 piOpt = cell(M, 1);
1761 for k = 1:M
1762 piOpt{k} = sol(piIdx{k})';
1763 end
1764else
1765 xOpt = [];
1766 piOpt = {};
1767end
1768
1769end
1770
1771%% Zero Potential Relaxation (zpr)
1772% Uses potential matrix to add cutting plane constraints
1773function [xOpt, piOpt, fval, exitflag] = solve_zpr(targetAction, boundType, ...
1774 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL)
1775
1776% First get initial solution using inap or from midpoint
1777x0 = (xL + xU) / 2;
1778[pi0, Q0] = compute_eq_zpr(x0, Aa, Pb, L, ACT, PSV, M, A, N);
1779
1780% Compute potential matrices using Jacobi iteration
1781g = cell(M, max(N));
1782for k = 1:M
1783 for n = 1:N(k)
1784 f_vec = zeros(N(k), 1);
1785 f_vec(n) = 1;
1786 g{k, n} = potential_jacobi(Q0{k}, pi0{k}, f_vec);
1787 end
1788end
1789
1790% Variables: [x(1:A); pi{1}(:); pi{2}(:); ...; z{1,1}(:); z{1,2}(:); ...]
1791nX = A;
1792nPi = sum(N);
1793nZ = 0;
1794for a = 1:A
1795 nZ = nZ + N(ACT(a)) + N(PSV(a));
1796end
1797nVar = nX + nPi + nZ;
1798
1799% Variable indexing
1800xIdx = 1:A;
1801piIdx = cell(M, 1);
1802offset = A;
1803for k = 1:M
1804 piIdx{k} = offset + (1:N(k));
1805 offset = offset + N(k);
1806end
1807zIdx = cell(A, 2);
1808for a = 1:A
1809 zIdx{a, 1} = offset + (1:N(ACT(a)));
1810 offset = offset + N(ACT(a));
1811 zIdx{a, 2} = offset + (1:N(PSV(a)));
1812 offset = offset + N(PSV(a));
1813end
1814
1815% Objective
1816f = zeros(nVar, 1);
1817if boundType == 0
1818 f(targetAction) = 1;
1819else
1820 f(targetAction) = -1;
1821end
1822
1823% Build constraints
1824Aeq = [];
1825beq = [];
1826Aineq = [];
1827bineq = [];
1828
1829% pi{k} * 1 = 1
1830for k = 1:M
1831 row = zeros(1, nVar);
1832 row(piIdx{k}) = 1;
1833 Aeq = [Aeq; row];
1834 beq = [beq; 1];
1835end
1836
1837% sum(z{a,type}) = x(a)
1838for a = 1:A
1839 row = zeros(1, nVar);
1840 row(zIdx{a, 1}) = 1;
1841 row(xIdx(a)) = -1;
1842 Aeq = [Aeq; row];
1843 beq = [beq; 0];
1844
1845 row = zeros(1, nVar);
1846 row(zIdx{a, 2}) = 1;
1847 row(xIdx(a)) = -1;
1848 Aeq = [Aeq; row];
1849 beq = [beq; 0];
1850end
1851
1852% pi{k} * Q{k} = 0 (balance equations)
1853for k = 1:M
1854 Qlocal = L{k} - diag(sum(L{k}, 2));
1855
1856 for s = 1:N(k)
1857 row = zeros(1, nVar);
1858 row(piIdx{k}) = Qlocal(:, s)';
1859
1860 for a = 1:A
1861 if PSV(a) == k
1862 Qp = Pb{a} - diag(sum(Pb{a}, 2));
1863 row(zIdx{a, 2}) = row(zIdx{a, 2}) + Qp(:, s)';
1864 elseif ACT(a) == k
1865 Qa = Aa{a} - diag(sum(Aa{a}, 2));
1866 row(piIdx{k}) = row(piIdx{k}) + Qa(:, s)';
1867 end
1868 end
1869
1870 Aineq = [Aineq; row; -row];
1871 bineq = [bineq; PFTOL; PFTOL];
1872 end
1873end
1874
1875% RC3 constraint: z{a,active}*I - pi*Aa = 0
1876for a = 1:A
1877 k = ACT(a);
1878 for i = 1:N(k)
1879 row = zeros(1, nVar);
1880 row(zIdx{a, 1}(i)) = 1;
1881 row(piIdx{k}) = -Aa{a}(:, i)';
1882
1883 Aineq = [Aineq; row; -row];
1884 bineq = [bineq; PFTOL; PFTOL];
1885 end
1886end
1887
1888% Potential constraints (cutting planes from zpr)
1889for k = 1:M
1890 for n = 1:N(k)
1891 % delta constraint: pi0 - pi + correction terms = 0
1892 row = zeros(1, nVar);
1893 row(piIdx{k}(n)) = -1;
1894 rhs = -pi0{k}(n);
1895
1896 for b = 1:A
1897 if PSV(b) == k
1898 Qp = Pb{b} - diag(sum(Pb{b}, 2));
1899 % z{b,passive} * Qp * g{k,n}
1900 coeff = Qp * g{k, n};
1901 row(zIdx{b, 2}) = row(zIdx{b, 2}) + coeff';
1902 % -x0(b) * pi * Qp * g{k,n}
1903 coeff2 = Qp * g{k, n};
1904 row(piIdx{k}) = row(piIdx{k}) - x0(b) * coeff2';
1905 end
1906 end
1907
1908 Aineq = [Aineq; row; -row];
1909 bineq = [bineq; rhs + PFTOL; -rhs + PFTOL];
1910 end
1911end
1912
1913% McCormick envelopes
1914for a = 1:A
1915 k_act = ACT(a);
1916 k_psv = PSV(a);
1917
1918 for i = 1:N(k_act)
1919 xLa = xL(a); xUa = xU(a);
1920 pL = piL{k_act}(i); pU = piU{k_act}(i);
1921
1922 row = zeros(1, nVar);
1923 row(zIdx{a, 1}(i)) = -1;
1924 row(xIdx(a)) = pL;
1925 row(piIdx{k_act}(i)) = xLa;
1926 Aineq = [Aineq; row];
1927 bineq = [bineq; pL * xLa];
1928
1929 row = zeros(1, nVar);
1930 row(zIdx{a, 1}(i)) = 1;
1931 row(xIdx(a)) = -pL;
1932 row(piIdx{k_act}(i)) = -xUa;
1933 Aineq = [Aineq; row];
1934 bineq = [bineq; -pL * xUa];
1935
1936 row = zeros(1, nVar);
1937 row(zIdx{a, 1}(i)) = 1;
1938 row(xIdx(a)) = -pU;
1939 row(piIdx{k_act}(i)) = -xLa;
1940 Aineq = [Aineq; row];
1941 bineq = [bineq; -pU * xLa];
1942
1943 row = zeros(1, nVar);
1944 row(zIdx{a, 1}(i)) = -1;
1945 row(xIdx(a)) = pU;
1946 row(piIdx{k_act}(i)) = xUa;
1947 Aineq = [Aineq; row];
1948 bineq = [bineq; pU * xUa];
1949 end
1950
1951 for i = 1:N(k_psv)
1952 xLa = xL(a); xUa = xU(a);
1953 pL = piL{k_psv}(i); pU = piU{k_psv}(i);
1954
1955 row = zeros(1, nVar);
1956 row(zIdx{a, 2}(i)) = -1;
1957 row(xIdx(a)) = pL;
1958 row(piIdx{k_psv}(i)) = xLa;
1959 Aineq = [Aineq; row];
1960 bineq = [bineq; pL * xLa];
1961
1962 row = zeros(1, nVar);
1963 row(zIdx{a, 2}(i)) = 1;
1964 row(xIdx(a)) = -pL;
1965 row(piIdx{k_psv}(i)) = -xUa;
1966 Aineq = [Aineq; row];
1967 bineq = [bineq; -pL * xUa];
1968
1969 row = zeros(1, nVar);
1970 row(zIdx{a, 2}(i)) = 1;
1971 row(xIdx(a)) = -pU;
1972 row(piIdx{k_psv}(i)) = -xLa;
1973 Aineq = [Aineq; row];
1974 bineq = [bineq; -pU * xLa];
1975
1976 row = zeros(1, nVar);
1977 row(zIdx{a, 2}(i)) = -1;
1978 row(xIdx(a)) = pU;
1979 row(piIdx{k_psv}(i)) = xUa;
1980 Aineq = [Aineq; row];
1981 bineq = [bineq; pU * xUa];
1982 end
1983end
1984
1985% Bounds
1986lb = zeros(nVar, 1);
1987ub = Inf(nVar, 1);
1988lb(xIdx) = xL;
1989ub(xIdx) = xU;
1990for k = 1:M
1991 lb(piIdx{k}) = piL{k};
1992 ub(piIdx{k}) = piU{k};
1993end
1994for a = 1:A
1995 lb(zIdx{a, 1}) = 0;
1996 ub(zIdx{a, 1}) = xU(a);
1997 lb(zIdx{a, 2}) = 0;
1998 ub(zIdx{a, 2}) = xU(a);
1999end
2000
2001% Solve LP
2002options = optimoptions('linprog', 'Algorithm', 'interior-point', 'Display', 'off');
2003[sol, fval, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, options);
2004
2005if exitflag > 0
2006 xOpt = sol(xIdx);
2007 piOpt = cell(M, 1);
2008 for k = 1:M
2009 piOpt{k} = sol(piIdx{k})';
2010 end
2011else
2012 xOpt = [];
2013 piOpt = {};
2014end
2015
2016end
2017
2018%% Tightened Zero Potential Relaxation (tzpr0inf/tzpr1inf)
2019% Combines zero potential relaxation with tightened LP constraints
2020% U = tightening level (0 for tzpr0inf, 1 for tzpr1inf)
2021function [xOpt, piOpt, fval, exitflag] = solve_tzpr(targetAction, boundType, ...
2022 Aa, Pb, L, ACT, PSV, M, A, N, xL, xU, piL, piU, PFTOL, U)
2023
2024% First get initial solution from midpoint
2025x0 = (xL + xU) / 2;
2026[pi0, Q0] = compute_eq_zpr(x0, Aa, Pb, L, ACT, PSV, M, A, N);
2027
2028% Compute potential matrices
2029g = cell(M, max(N));
2030for k = 1:M
2031 for n = 1:N(k)
2032 f_vec = zeros(N(k), 1);
2033 f_vec(n) = 1;
2034 g{k, n} = potential_jacobi(Q0{k}, pi0{k}, f_vec);
2035 end
2036end
2037
2038% Variables: [x(1:A); pi{1}(:); pi{2}(:); ...; z{1,1}(:); z{1,2}(:); ...]
2039nX = A;
2040nPi = sum(N);
2041nZ = 0;
2042for a = 1:A
2043 nZ = nZ + N(ACT(a)) + N(PSV(a));
2044end
2045nVar = nX + nPi + nZ;
2046
2047% Variable indexing
2048xIdx = 1:A;
2049piIdx = cell(M, 1);
2050offset = A;
2051for k = 1:M
2052 piIdx{k} = offset + (1:N(k));
2053 offset = offset + N(k);
2054end
2055zIdx = cell(A, 2);
2056for a = 1:A
2057 zIdx{a, 1} = offset + (1:N(ACT(a)));
2058 offset = offset + N(ACT(a));
2059 zIdx{a, 2} = offset + (1:N(PSV(a)));
2060 offset = offset + N(PSV(a));
2061end
2062
2063% Objective
2064f = zeros(nVar, 1);
2065if boundType == 0
2066 f(targetAction) = 1;
2067else
2068 f(targetAction) = -1;
2069end
2070
2071% Build constraints
2072Aeq = [];
2073beq = [];
2074Aineq = [];
2075bineq = [];
2076
2077% pi{k} * 1 = 1
2078for k = 1:M
2079 row = zeros(1, nVar);
2080 row(piIdx{k}) = 1;
2081 Aeq = [Aeq; row];
2082 beq = [beq; 1];
2083end
2084
2085% sum(z{a,type}) = x(a)
2086for a = 1:A
2087 row = zeros(1, nVar);
2088 row(zIdx{a, 1}) = 1;
2089 row(xIdx(a)) = -1;
2090 Aeq = [Aeq; row];
2091 beq = [beq; 0];
2092
2093 row = zeros(1, nVar);
2094 row(zIdx{a, 2}) = 1;
2095 row(xIdx(a)) = -1;
2096 Aeq = [Aeq; row];
2097 beq = [beq; 0];
2098end
2099
2100% pi{k} * Q{k} = 0 (balance equations)
2101for k = 1:M
2102 Qlocal = L{k} - diag(sum(L{k}, 2));
2103
2104 for s = 1:N(k)
2105 row = zeros(1, nVar);
2106 row(piIdx{k}) = Qlocal(:, s)';
2107
2108 for a = 1:A
2109 if PSV(a) == k
2110 Qp = Pb{a} - diag(sum(Pb{a}, 2));
2111 row(zIdx{a, 2}) = row(zIdx{a, 2}) + Qp(:, s)';
2112 elseif ACT(a) == k
2113 Qa = Aa{a} - diag(sum(Aa{a}, 2));
2114 row(piIdx{k}) = row(piIdx{k}) + Qa(:, s)';
2115 end
2116 end
2117
2118 Aineq = [Aineq; row; -row];
2119 bineq = [bineq; PFTOL; PFTOL];
2120 end
2121end
2122
2123% RC3 constraint
2124for a = 1:A
2125 k = ACT(a);
2126 for i = 1:N(k)
2127 row = zeros(1, nVar);
2128 row(zIdx{a, 1}(i)) = 1;
2129 row(piIdx{k}) = -Aa{a}(:, i)';
2130
2131 Aineq = [Aineq; row; -row];
2132 bineq = [bineq; PFTOL; PFTOL];
2133 end
2134end
2135
2136% Potential constraints (from zpr)
2137for k = 1:M
2138 for n = 1:N(k)
2139 row = zeros(1, nVar);
2140 row(piIdx{k}(n)) = -1;
2141 rhs = -pi0{k}(n);
2142
2143 for b = 1:A
2144 if PSV(b) == k
2145 Qp = Pb{b} - diag(sum(Pb{b}, 2));
2146 coeff = Qp * g{k, n};
2147 row(zIdx{b, 2}) = row(zIdx{b, 2}) + coeff';
2148 coeff2 = Qp * g{k, n};
2149 row(piIdx{k}) = row(piIdx{k}) - x0(b) * coeff2';
2150 end
2151 end
2152
2153 Aineq = [Aineq; row; -row];
2154 bineq = [bineq; rhs + PFTOL; -rhs + PFTOL];
2155 end
2156end
2157
2158% Tightened LP constraints (from tlpr)
2159% U is passed as parameter: 0 for tzpr0inf, 1 for tzpr1inf
2160for k = 1:M
2161 for b = 1:A
2162 if ACT(b) == k
2163 Qlocal = L{k} - diag(sum(L{k}, 2));
2164
2165 for u = 1:U
2166 Aa_pow = Aa{b}^u;
2167
2168 for s = 1:N(k)
2169 row = zeros(1, nVar);
2170 coeffZ = (Aa{b}^(u-1)) * Qlocal;
2171 row(zIdx{b, 1}) = coeffZ(:, s)';
2172
2173 for c = 1:A
2174 if PSV(c) == k
2175 Qp = Pb{c} - diag(sum(Pb{c}, 2));
2176 coeffZp = Aa_pow * Qp;
2177 row(zIdx{c, 2}) = row(zIdx{c, 2}) + coeffZp(:, s)';
2178 elseif ACT(c) == k
2179 Qa = Aa{c} - diag(sum(Aa{c}, 2));
2180 coeffZa = (Aa{b}^(u-1)) * Qa;
2181 row(zIdx{b, 1}) = row(zIdx{b, 1}) + coeffZa(:, s)';
2182 end
2183 end
2184
2185 Aineq = [Aineq; row; -row];
2186 bineq = [bineq; PFTOL; PFTOL];
2187 end
2188 end
2189
2190 % Infinity constraint
2191 try
2192 Aa_inv = inv(eye(N(k)) - Aa{b});
2193
2194 for s = 1:N(k)
2195 row = zeros(1, nVar);
2196 coeffZ = Aa_inv * Qlocal;
2197 row(zIdx{b, 1}) = coeffZ(:, s)';
2198
2199 for c = 1:A
2200 if PSV(c) == k
2201 Qp = Pb{c} - diag(sum(Pb{c}, 2));
2202 coeffZp = Aa{b} * Aa_inv * Qp;
2203 row(zIdx{c, 2}) = row(zIdx{c, 2}) + coeffZp(:, s)';
2204 elseif ACT(c) == k
2205 Qa = Aa{c} - diag(sum(Aa{c}, 2));
2206 coeffZa = Aa_inv * Qa;
2207 row(zIdx{b, 1}) = row(zIdx{b, 1}) + coeffZa(:, s)';
2208 end
2209 end
2210
2211 Aineq = [Aineq; row; -row];
2212 bineq = [bineq; PFTOL; PFTOL];
2213 end
2214 catch
2215 % Matrix not invertible
2216 end
2217 end
2218 end
2219end
2220
2221% McCormick envelopes
2222for a = 1:A
2223 k_act = ACT(a);
2224 k_psv = PSV(a);
2225
2226 for i = 1:N(k_act)
2227 xLa = xL(a); xUa = xU(a);
2228 pL = piL{k_act}(i); pU = piU{k_act}(i);
2229
2230 row = zeros(1, nVar);
2231 row(zIdx{a, 1}(i)) = -1;
2232 row(xIdx(a)) = pL;
2233 row(piIdx{k_act}(i)) = xLa;
2234 Aineq = [Aineq; row];
2235 bineq = [bineq; pL * xLa];
2236
2237 row = zeros(1, nVar);
2238 row(zIdx{a, 1}(i)) = 1;
2239 row(xIdx(a)) = -pL;
2240 row(piIdx{k_act}(i)) = -xUa;
2241 Aineq = [Aineq; row];
2242 bineq = [bineq; -pL * xUa];
2243
2244 row = zeros(1, nVar);
2245 row(zIdx{a, 1}(i)) = 1;
2246 row(xIdx(a)) = -pU;
2247 row(piIdx{k_act}(i)) = -xLa;
2248 Aineq = [Aineq; row];
2249 bineq = [bineq; -pU * xLa];
2250
2251 row = zeros(1, nVar);
2252 row(zIdx{a, 1}(i)) = -1;
2253 row(xIdx(a)) = pU;
2254 row(piIdx{k_act}(i)) = xUa;
2255 Aineq = [Aineq; row];
2256 bineq = [bineq; pU * xUa];
2257 end
2258
2259 for i = 1:N(k_psv)
2260 xLa = xL(a); xUa = xU(a);
2261 pL = piL{k_psv}(i); pU = piU{k_psv}(i);
2262
2263 row = zeros(1, nVar);
2264 row(zIdx{a, 2}(i)) = -1;
2265 row(xIdx(a)) = pL;
2266 row(piIdx{k_psv}(i)) = xLa;
2267 Aineq = [Aineq; row];
2268 bineq = [bineq; pL * xLa];
2269
2270 row = zeros(1, nVar);
2271 row(zIdx{a, 2}(i)) = 1;
2272 row(xIdx(a)) = -pL;
2273 row(piIdx{k_psv}(i)) = -xUa;
2274 Aineq = [Aineq; row];
2275 bineq = [bineq; -pL * xUa];
2276
2277 row = zeros(1, nVar);
2278 row(zIdx{a, 2}(i)) = 1;
2279 row(xIdx(a)) = -pU;
2280 row(piIdx{k_psv}(i)) = -xLa;
2281 Aineq = [Aineq; row];
2282 bineq = [bineq; -pU * xLa];
2283
2284 row = zeros(1, nVar);
2285 row(zIdx{a, 2}(i)) = -1;
2286 row(xIdx(a)) = pU;
2287 row(piIdx{k_psv}(i)) = xUa;
2288 Aineq = [Aineq; row];
2289 bineq = [bineq; pU * xUa];
2290 end
2291end
2292
2293% Bounds
2294lb = zeros(nVar, 1);
2295ub = Inf(nVar, 1);
2296lb(xIdx) = xL;
2297ub(xIdx) = xU;
2298for k = 1:M
2299 lb(piIdx{k}) = piL{k};
2300 ub(piIdx{k}) = piU{k};
2301end
2302for a = 1:A
2303 lb(zIdx{a, 1}) = 0;
2304 ub(zIdx{a, 1}) = xU(a);
2305 lb(zIdx{a, 2}) = 0;
2306 ub(zIdx{a, 2}) = xU(a);
2307end
2308
2309% Solve LP
2310options = optimoptions('linprog', 'Algorithm', 'interior-point', 'Display', 'off');
2311[sol, fval, exitflag] = linprog(f, Aineq, bineq, Aeq, beq, lb, ub, options);
2312
2313if exitflag > 0
2314 xOpt = sol(xIdx);
2315 piOpt = cell(M, 1);
2316 for k = 1:M
2317 piOpt{k} = sol(piIdx{k})';
2318 end
2319else
2320 xOpt = [];
2321 piOpt = {};
2322end
2323
2324end
2325
2326%% Compute equilibrium for zpr
2327function [pi, Q] = compute_eq_zpr(x, Aa, Pb, L, ACT, PSV, M, A, N)
2328Q = cell(1, M);
2329pi = cell(1, M);
2330
2331for k = 1:M
2332 Qk = L{k} - diag(sum(L{k}, 2));
2333
2334 for a = 1:A
2335 if PSV(a) == k
2336 Qk = Qk + x(a) * (Pb{a} - diag(sum(Pb{a}, 2)));
2337 elseif ACT(a) == k
2338 Qk = Qk + Aa{a} - diag(sum(Aa{a}, 2));
2339 end
2340 end
2341
2342 Q{k} = Qk;
2343 Qk_inf = ctmc_makeinfgen(Qk);
2344 pi{k} = ctmc_solve(Qk_inf);
2345end
2346end
2347
2348%% Potential matrix computation using Jacobi iteration
2349function g = potential_jacobi(Q, pi, f)
2350% Compute potential g such that Q*g = f - pi*f*e
2351% Uses Jacobi iteration method
2352
2353n = length(pi);
2354tol = 1e-10;
2355maxIter = 1000;
2356
2357% Right-hand side
2358b = f - (pi * f) * ones(n, 1);
2359
2360% Initialize g
2361g = zeros(n, 1);
2362
2363% Extract diagonal
2364d = diag(Q);
2365
2366% Jacobi iteration
2367for iter = 1:maxIter
2368 g_old = g;
2369
2370 for i = 1:n
2371 if abs(d(i)) > 1e-12
2372 sum_term = Q(i, :) * g_old - d(i) * g_old(i);
2373 g(i) = (b(i) - sum_term) / d(i);
2374 end
2375 end
2376
2377 if norm(g - g_old, inf) < tol
2378 break;
2379 end
2380end
2381
2382% Normalize to satisfy pi*g = 0
2383g = g - (pi * g);
2384
2385end
Generated by doxygen 1.9.8