Garicano(2000)模型符号优化求解：Matlab/Python实现遇阻求助

我正尝试复现Garicano(2000)中的优化问题，目标是最大化指定方程，求解$z_w$和$z_m$以得到目标解析解。我分别尝试了Matlab和Python实现，但均遇到问题：

Matlab实现问题 #

代码如下：

%%   Model of firm layers' with knowledge transmission & Trainment
clear all;
clc;

syms z_w z_m z h t w n a1 a2 a3 L

% Here we are using the objective function with the restriction included

F = ((1-exp(-L*z_m))/(h*(1-(1-exp(-L*z_w))))) - (t - a3)*z_m - w*(a3) - (((t - a1).*z_w + w*(a1))/(h*(1-exp(-L*z_w))));


solve(F,z_w)

运行结果提示：

Warning: Unable to find explicit solution. For options, see help. 
In sym/solve (line 317)
In untitled3 (line 12) 
 
ans = 

Empty sym: 0-by-1

我推测这是因无数值输入导致，但不确定解决方法。

Python实现问题 #

代码如下：

from __future__ import division
from simpy import *
from sympy import symbols, exp, diff, solve

### We define the variables and parameters

z_m, z_w, z = symbols('z_m z_w z', positive = True)
L = symbols('lambda')
h, t, w, n, a1, a2, a3, L = symbols('h t w n a1 a2 a3 L', Real = True)
F, C1, m = symbols('F C1 m')

# Defining the function
Objective = ((1-exp(-L*z_m))/(h*(1-(1-exp(-L*z_w)))))-(t - a3)*z_m-w*(a3)-(((t - a1)*z_w + w*(a1))/(h*(1-exp(-L*z_w))))
Objective

C1 = h*n*(1-exp(m-L*z_w))
C1

# Identify the Objective function
Lag = Objective
Lag

#First order conditions

FOC_1 = diff(Lag,z_w)
FOC_2 = diff(Lag,z_m)


Solutions = solve([FOC_1, FOC_2], [z_w, z_m])

运行后抛出错误：

---------------------------------------------------------------------------
NotImplementedError                       Traceback (most recent call last)
<ipython-input-8-9cff7d4c192f> in <module>
----> 1 Solutions = solve([FOC_1, FOC_2], [z_w, z_m])

~\anaconda3\lib\site-packages\sympy\solvers\solvers.py in solve(f, *symbols, **flags)
   1094         solution = _solve(f[0], *symbols, **flags)
   1095     else:
-> 1096         solution = _solve_system(f, symbols, **flags)
   1097 
   1098     #

~\anaconda3\lib\site-packages\sympy\solvers\solvers.py in _solve_system(exprs, symbols, **flags)
   1887                     got_s.add(s)
   1888                 if not hit:
-> 1889                     raise NotImplementedError('could not solve %s' % eq2)
   1890             else:
   1891                 result = newresult

NotImplementedError: could not solve h**4*(-1 + h*(-a3 + t)*exp(L*z_m)/L)**3*(-a3 + t)**2*(exp(L*z_m) - 1)*exp(2*L*z_m)/L + h**4*(-1 + h*(-a3 + t)*exp(L*z_m)/L)*(-a3 + t)**2*(a1*w + (-a1 + t)*(L*z_m + log(h*(-a3 + t)/L))/L)*exp(3*L*z_m)/L - h**4*(-1 + h*(-a3 + t)*exp(L*z_m)/L)**2*(-a1 + t)*(-a3 + t)**2*exp(3*L*z_m)/L**2

我对正确构建符号优化问题理解不足，恳请提供帮助，任何反馈都将十分感谢。我已尝试调整代码使用符号函数，但仍无法解决问题。

内容的提问来源于stack exchange，提问作者Daniel Garzon Hernandez