在 Matlab 中最大化函数（使用 fsolve 和导数）

Question

我想找到向量p（即p(1),p(2),p(3),...）的值，它可以最大化函数A(p)。

我正在使用 MATLAB 来执行此操作，我发现 fsolve 我认为它可以帮助我。所以我做了函数A：

function A = myfun(p)

    R  = 0.1; 
    u1 = 500;
    u2 = 400;
    u3 = 300;

    A = ( (p(1)+p(2)+p(3)) * (1/u1+1/u2+1/u3)) * ...
        (1 + R*(p(1)^2+p(2)^2+p(3)^2) * (1/u1+1/u2+1/u3) );

然后我需要求解一个方程组：

diff(A,p(1))==0

diff(A,p(2))==0

diff(A,p(3))==0

生成的p 向量将解决我的问题。

fsolve 如何解决这个方程组 (p0=[1 1 1])？

Answer 1

下面是一个使用符号工具箱的例子：

%// (there's probably a way to generalize this as well)
dAdp = cellstr(char(...
    [char(diff('C * (p1+p2+p3) * (1 + R*C*(p1^2+p2^2+p3^2))', 'p1')) ';'],...
    [char(diff('C * (p1+p2+p3) * (1 + R*C*(p1^2+p2^2+p3^2))', 'p2')) ';'],...
    [char(diff('C * (p1+p2+p3) * (1 + R*C*(p1^2+p2^2+p3^2))', 'p3')) ';']))

%// convert to proper vector equation
dAdp = regexprep(dAdp, 'p([0-9])', 'p\($1\)');

%// convert to function handle
F = str2func( strcat('@(p) [', dAdp{:}, ']') );

但我不建议这样做（它完全不可读且容易出错）。

您可以编写一个适当的函数并使用符号工具箱评估上面的dAdp，但我也不建议这样做（它非常慢）。

我是一个数字人，因为我只相信计算机的用途：计算。我在纸上进行推导，除了过于简单但冗长乏味的推导（您可能会说这是这样的情况，我仍然喜欢在纸上进行，因为我也喜欢锻炼我的脑 ：）。

我建议您这样做，和/或使用符号工具箱不断检查自己。恕我直言，您应该将其用作辅助，而不是作为主要引擎。

那么，我们开始吧。这又是你的函数，这次以不同的形式，应用了更多的大脑：

function A = myfun(p)    
    R = 0.1; 
    u = [500; 400; 300];     
    C = sum(1./u);
    B = C * sum(p) * (1 + R*C*sum(p.^2));

所以，你想解决 ∇A(p) = 0。最好的办法是多动脑筋。您可以验证向量导数等于：

function F = Aprime(p)
    R = 0.1;
    u = [500; 400; 300];
    C = sum(1./u);
    F = C*( C*R*sum(p.^2) + 2*C*R*p*sum(p) + 1 );

你可以用更多的大脑来解决：

C·( CR·Σp² + 2CR·p·Σp + 1 ) = 0
                Σp² + 2p·Σp = -1/(CR)

vector == scalar：这意味着 p 中的所有元素都是相等的。
代入q = p1 = p2 = p3，则

 3q² + 2q·3q = -1/(CR)
         9q² = -1/(CR)

⇒ q = ⅓·√(-1/(CR))

（表示有问题），或者在fsolve 中这样：

fsolve(@Aprime, [1 1 1])

这将立即导致

fsolve 完成，因为函数值的向量在初始 点接近于零，由函数的默认值测量 容差，并且问题看起来是有规律的 渐变。

这确实表示麻烦。

既然你已经表明你对最大值感兴趣，并且函数似乎没有固定点，唯一的结论是函数没有最大值也没有最小值，除非你对 p 施加界限。事实上，如果你将维数减少到 2，并绘制一个图：

R = 0.1;
u = [500; 400; 300];
C = sum(1./u);
B = @(p) C * sum(p) .* (1 + R*C*sum(p.^2));

[p1,p2] = meshgrid(-10:0.1:10);
surf(p1,p2,reshape(B([p1(:) p2(:)].'), size(p1)), 'edgecolor', 'none')

...没有极值。