Iteration method of optimization

I am trying to solve the problem using stackelberg sequential equilibrium method -

Understanding the Approach:

In a Stackelberg game, we have two players: the leader and the follower. The leader chooses its strategy first, knowing that the follower will then optimize their strategy in response to the leader's choice. This sequential decision-making process requires us to iteratively solve for the best responses of each player until we converge to a Nash equilibrium.

-Follower's Best Response (t given k):

The follower's strategy t is a function of the leader's strategy k. We need to find the value of t that maximizes the follower's utility given the leader's choice of k. We start by initial guess of k and compute corresponding optimal t.

Leader's Best Response (k inputing response of t):

Once we have the optimal t, the leader then chooses its strategy k to maximize its utility, knowing the follower will choose t optimally in response.

We will do this iteration till initial guess converge to optimal k. Following code encapsulate my idea -

% Parameters
m_values = linspace(0, 1, 100); % Range of m values
tolerance = 1e-6; % Convergence tolerance
max_iter = 1000; % Maximum number of iterations
 % Initialize arrays to store results
k_solutions = zeros(size(m_values));
t_solutions = zeros(size(m_values));
 % Function to compute follower's best response t given k and m
compute_t = @(k, m) fminbnd(@(t) abs(sqrt(k) * ((m./t + 2) / 3) - sqrt(1 - k) * (1 / 3) * (2 + m./t + (2 * m + t) ./ ((m - t).^2 - 2 * t))), 0, 1);
 % Function to compute leader's objective function given k, t, and m
compute_obj = @(k, t, m) -(sqrt(1 - k) * (t^6 * (2 * m * (3 * m + 2) + 4) - t^5 * (m * (m * (5 * m + 6) + 11) + 12) + m^6 + 2 * t^7 - t^8 - t^4 * (m * (3 * m * (m * (3 * m + 2) - 3) - 16) - 9) - 2 * m^3 * t^2 * (3 * m^3 + 2 * m + 3) + m^2 * t^3 * (m * (3 * m * (5 * m + 4) + 1) - 4)) + sqrt(k) * (t^4 * (m * (9 * m * (m * (m + 4) + 5) + 28) + 9) - m^4 * (3 * m - 2) - t^6 * (6 * m * (m + 2) + 10) + t^5 * (4 * m^3 + 12 * m + 8) + t^8 + m * t^2 * (m * (m * (4 * m^3 + 3 * m + 9) + 6) - 8) + 2 * m^2 * t * (3 * m - 2) - 2 * m * t^3 * (m * (m * (6 * m * (m + 2) + 7) + 6) - 2))) / (18 * t^3 * (2 * t - (m - t)^2));
 % Main loop to solve for each m
for i = 1:length(m_values)
 m = m_values(i);
 k_guess = 0.75; % Initial guess for k
 k_opt = k_guess;
 t_opt = compute_t(k_opt, m);
 
 % Iterative optimization
 for iter = 1:max_iter
 % Compute optimal t given k
 t_opt = compute_t(k_opt, m);
 
 % Optimize k given t
 k_opt_new = fminbnd(@(k) compute_obj(k, t_opt, m), 0.5, 1); % Use fminbnd for bounded optimization
 
 % Check for convergence
 if abs(k_opt_new - k_opt) < tolerance
 k_opt = k_opt_new;
 break;
 end
 k_opt = k_opt_new;
 end
 
 % Store solutions
 k_solutions(i) = k_opt;
 t_solutions(i) = t_opt;
end
 % Plot results
figure;
plot(m_values, k_solutions, 'b-', 'LineWidth', 1.5);
hold on;
plot(m_values, t_solutions, 'r--', 'LineWidth', 1.5);
xlabel('m');
ylabel('Value');
legend('Optimal k', 'Optimal t');
title('Stackelberg Equilibrium Solutions');

 % Display final values
fprintf('Final solutions:\n');

Final solutions:

fprintf('m \t k \t t\n');

m 	 k 	 t

for i = 1:length(m_values)
 fprintf('%.4f \t %.4f \t %.4f\n', m_values(i), k_solutions(i), t_solutions(i));
end

0.0000 	 0.5001 	 0.0001
0.0101 	 0.6188 	 0.2044
0.0202 	 0.7234 	 0.3086
0.0303 	 0.7870 	 0.3942
0.0404 	 0.8288 	 0.4686
0.0505 	 0.8575 	 0.5345
0.0606 	 0.8775 	 0.5934
0.0707 	 0.8919 	 0.6464
0.0808 	 0.9022 	 0.6940
0.0909 	 0.9099 	 0.7372
0.1010 	 0.9155 	 0.7766
0.1111 	 0.9198 	 0.8126
0.1212 	 0.9229 	 0.8458
0.1313 	 0.9253 	 0.8767
0.1414 	 0.9271 	 0.9054
0.1515 	 0.9284 	 0.9324
0.1616 	 0.9293 	 0.9578
0.1717 	 0.9300 	 0.9821
0.1818 	 0.9304 	 1.0000
0.1919 	 0.9309 	 1.0000
0.2020 	 0.9314 	 1.0000
0.2121 	 0.9321 	 1.0000
0.2222 	 0.9328 	 1.0000
0.2323 	 0.9335 	 1.0000
0.2424 	 0.9344 	 1.0000
0.2525 	 0.9353 	 1.0000
0.2626 	 0.9363 	 1.0000
0.2727 	 0.9373 	 1.0000
0.2828 	 0.9383 	 1.0000
0.2929 	 0.9394 	 1.0000
0.3030 	 0.9405 	 1.0000
0.3131 	 0.9417 	 1.0000
0.3232 	 0.9428 	 1.0000
0.3333 	 0.9441 	 1.0000
0.3434 	 0.9453 	 1.0000
0.3535 	 0.9465 	 0.9999
0.3636 	 0.9477 	 0.9999
0.3737 	 0.9490 	 0.9999
0.3838 	 0.9502 	 0.9999
0.3939 	 0.9514 	 0.9999
0.4040 	 0.9527 	 0.9999
0.4141 	 0.9539 	 0.9999
0.4242 	 0.9551 	 0.9999
0.4343 	 0.9563 	 0.9999
0.4444 	 0.9575 	 0.9999
0.4545 	 0.9587 	 0.9999
0.4646 	 0.9599 	 0.9999
0.4747 	 0.9610 	 0.9999
0.4848 	 0.9621 	 0.9999
0.4949 	 0.9632 	 0.9999
0.5051 	 0.9643 	 0.9999
0.5152 	 0.9653 	 0.9999
0.5253 	 0.9664 	 0.9999
0.5354 	 0.9674 	 0.9999
0.5455 	 0.9683 	 0.9999
0.5556 	 0.9693 	 0.9999
0.5657 	 0.9702 	 0.9999
0.5758 	 0.9711 	 0.9999
0.5859 	 0.9720 	 0.9999
0.5960 	 0.9728 	 0.9999
0.6061 	 0.9736 	 0.9999
0.6162 	 0.9744 	 0.9999
0.6263 	 0.9752 	 0.9999
0.6364 	 0.9759 	 0.9999
0.6465 	 0.9766 	 0.9999
0.6566 	 0.9773 	 0.9999
0.6667 	 0.9780 	 0.9999
0.6768 	 0.9786 	 0.9999
0.6869 	 0.9792 	 0.9999
0.6970 	 0.9798 	 0.9999
0.7071 	 0.9803 	 0.9999
0.7172 	 0.9809 	 0.9999
0.7273 	 0.9814 	 0.9999
0.7374 	 0.9819 	 0.9999
0.7475 	 0.9823 	 0.9999
0.7576 	 0.9828 	 0.9999
0.7677 	 0.9832 	 0.9999
0.7778 	 0.9836 	 0.9999
0.7879 	 0.9839 	 0.9999
0.7980 	 0.9843 	 0.9999
0.8081 	 0.9847 	 0.9999
0.8182 	 0.9850 	 1.0000
0.8283 	 0.9853 	 1.0000
0.8384 	 0.9856 	 1.0000
0.8485 	 0.9858 	 1.0000
0.8586 	 0.9861 	 1.0000
0.8687 	 0.9863 	 1.0000
0.8788 	 0.9865 	 1.0000
0.8889 	 0.9867 	 1.0000
0.8990 	 0.9869 	 1.0000
0.9091 	 0.9870 	 1.0000
0.9192 	 0.9872 	 1.0000
0.9293 	 0.9873 	 1.0000
0.9394 	 0.9874 	 1.0000
0.9495 	 0.9875 	 1.0000
0.9596 	 0.9876 	 1.0000
0.9697 	 0.9877 	 1.0000
0.9798 	 0.9878 	 1.0000
0.9899 	 0.9878 	 1.0000
1.0000 	 0.9878 	 1.0000

I could have used calculus... But since expression of t cannot be explicitly expressed as a function of k, substituting t into k and then differentiating wrt k (where, t is also function k) would give an untractably long expression.

Is the Approach I used correct? Can someone please suggest a better approach? Anyone please help

[Sabrina Garland]

Any help @Torsten

[Sabrina Garland]

Any help @anyone. Please suggest a way to improve the code

[Sabrina Garland]

@Torsten Please provide an insight

Answer 1

[Torsten]

That's what I get from your description. But I might have misunderstood something.

% Parameters
m_values = linspace(0, 1, 100); % Range of m values
% Function to compute leader's objective function given k, t, and m
compute_k_from_t = @(k, t, m) -(sqrt(1 - k) * (t^6 * (2 * m * (3 * m + 2) + 4) - t^5 * (m * (m * (5 * m + 6) + 11) + 12) + m^6 + 2 * t^7 - t^8 - t^4 * (m * (3 * m * (m * (3 * m + 2) - 3) - 16) - 9) - 2 * m^3 * t^2 * (3 * m^3 + 2 * m + 3) + m^2 * t^3 * (m * (3 * m * (5 * m + 4) + 1) - 4)) + sqrt(k) * (t^4 * (m * (9 * m * (m * (m + 4) + 5) + 28) + 9) - m^4 * (3 * m - 2) - t^6 * (6 * m * (m + 2) + 10) + t^5 * (4 * m^3 + 12 * m + 8) + t^8 + m * t^2 * (m * (m * (4 * m^3 + 3 * m + 9) + 6) - 8) + 2 * m^2 * t * (3 * m - 2) - 2 * m * t^3 * (m * (m * (6 * m * (m + 2) + 7) + 6) - 2))) / (18 * t^3 * (2 * t - (m - t)^2));
% Function to compute follower's best response t given k and m
compute_t_from_k = @(k, t, m) abs(sqrt(k) * ((m./t + 2) / 3) - sqrt(1 - k) * (1 / 3) * (2 + m./t + (2 * m + t) ./ ((m - t).^2 - 2 * t)));
% Main loop to solve for each m
k_values = linspace(0.5,1,1000);
for i = 1:numel(m_values)
  m = m_values(i);
  for j = 1:numel(k_values)
    k = k_values(j);
    t(j) = fminbnd(@(t)compute_t_from_k(k,t,m),0,1);
    leaders_gain(j) = compute_k_from_t(k,t(j),m);
  end
  [~,idx] = min(leaders_gain);
  k_opt(i) = k_values(idx);
  t_opt(i) = t(idx);
end
plot(m_values,[k_opt;t_opt])

[Sabrina Garland]

Thank you for your response @Torsten. The idea of looping over all values of k in the range is just brilliant! Instead of min(leader_gain), it should be maximised but the result remain the same (I wonder why?). Thanks again!

[Torsten]

The function "compute_k_from_t" has a minus sign at the beginning. Thus I thought you want to maximize (sqrt(1-k) * ... and you took -(sqrt(1-k) * ... to be minimized in which case you still had to take the index for which -(sqrt(1-k) * ... is minimized in the command

[~,idx] = min(leaders_gain);

because this is the index for which (sqrt(1-k) * ... is maximized.

Everything clear ? :-)

And the solution curves dramatically change in case you take the maximum:

% Parameters
m_values = linspace(0, 1, 100); % Range of m values
% Function to compute leader's objective function given k, t, and m
compute_k_from_t = @(k, t, m) -(sqrt(1 - k) * (t^6 * (2 * m * (3 * m + 2) + 4) - t^5 * (m * (m * (5 * m + 6) + 11) + 12) + m^6 + 2 * t^7 - t^8 - t^4 * (m * (3 * m * (m * (3 * m + 2) - 3) - 16) - 9) - 2 * m^3 * t^2 * (3 * m^3 + 2 * m + 3) + m^2 * t^3 * (m * (3 * m * (5 * m + 4) + 1) - 4)) + sqrt(k) * (t^4 * (m * (9 * m * (m * (m + 4) + 5) + 28) + 9) - m^4 * (3 * m - 2) - t^6 * (6 * m * (m + 2) + 10) + t^5 * (4 * m^3 + 12 * m + 8) + t^8 + m * t^2 * (m * (m * (4 * m^3 + 3 * m + 9) + 6) - 8) + 2 * m^2 * t * (3 * m - 2) - 2 * m * t^3 * (m * (m * (6 * m * (m + 2) + 7) + 6) - 2))) / (18 * t^3 * (2 * t - (m - t)^2));
% Function to compute follower's best response t given k and m
compute_t_from_k = @(k, t, m) abs(sqrt(k) * ((m./t + 2) / 3) - sqrt(1 - k) * (1 / 3) * (2 + m./t + (2 * m + t) ./ ((m - t).^2 - 2 * t)));
% Main loop to solve for each m
k_values = linspace(0.5,1,1000);
for i = 1:numel(m_values)
  m = m_values(i);
  for j = 1:numel(k_values)
    k = k_values(j);
    t(j) = fminbnd(@(t)compute_t_from_k(k,t,m),0,1);
    leaders_gain(j) = compute_k_from_t(k,t(j),m);
  end
  [~,idx] = max(leaders_gain);
  k_opt(i) = k_values(idx);
  t_opt(i) = t(idx);
end
plot(m_values,[k_opt;t_opt])

[Sabrina Garland]

Oops... Right... I got it... I took maximum after removing minus sign... That was a stupid question. Thank you for explaining it so well.

[Sabrina Garland]

Just last question @Torsten. Can something be wrong as we get solution to be more erratic than following a strict upward or downward path. Is there a way to check why solution are so erratic for lower values of m?

[Torsten]

I'm a bit surprised that the follower wants to minimize "compute_t_from_k". Is this a cost function and not a profit function ?

I'd suggest you test the method for a problem you know the solution of.

[Sabrina Garland]

t is kind of indisputable baggage... I don't have any running model to check my claims but let me try with tweaking things around. Thank you for all your help

Answer 2

[Umar]

Hi Sabrina, Based on the detailed explanation and code snippet provided, it seems that you are on the right track with your Stackelberg sequential equilibrium method approach. The approach you outlined involves iteratively solving for the best responses of the leader and follower until convergence to a Nash equilibrium is achieved. Your code snippet effectively captures the essence of the Stackelberg game setup, including the computation of the follower's best response given the leader's strategy and the leader's best response taking into account the follower's optimal strategy. The iterative optimization process you have implemented appears to be correctly structured to converge to optimal values of k and t for different values of m. While your current approach is valid, there are alternative methods that could potentially enhance efficiency or accuracy in solving for the Stackelberg equilibrium solutions. One suggestion could be to explore numerical optimization techniques beyond fminbnd, such as gradient-based methods like gradient descent or quasi-Newton methods like BFGS, which may offer faster convergence rates in some cases. Additionally, considering the complexity of the objective functions involved, you may also want to investigate symbolic computation libraries that could help in handling derivatives and optimizations more effectively. Libraries like SymPy in Python or symbolic math toolboxes in MATLAB can assist in dealing with complex mathematical expressions and enable more sophisticated analyses. In summary, while your current approach is sound, exploring alternative numerical optimization methods and symbolic computation tools could potentially provide insights into improving the efficiency and robustness of your solution method for solving Stackelberg games. Feel free to experiment with different techniques and tools to refine your approach further.