求助各位大佬，想绘制出上述函数表达式关于a的最大值的图像，其中想以b和s为横纵坐标，q作为动态调节的参数，绘制出这个函数表达式最大值高于某一个定值的b和s所确定的区域图像。目前只能对整个表达式求出关于a的最大值的带b、q、s参数的表达式，再给b、q、s赋值做regionplot画图，但是式子太复杂，带参数的最大值表达式计算不出来。请问各位大佬，有没有什么方法能够实时给三个参数赋值，而后再求出关于a的函数最大值的数值，再进行描点绘图呢。谢谢各位大佬！ 之前写的绘图代码如下
Manipulate[
RegionPlot[
MaxValue[{1/
2 a q (1 - 2 (-1 + a) b^2 + (
a^2 b s)/(-1 +
a) - ((2 + (-2 + a) a) s^2)/(-1 + a)^2 + ((-1 + a) b)/(
b - a b + s)),
0 < a < 1 && 0 < q && 0 < b < -(1/(-2 + a)) &&
0 < s < ((-1 + a) (1 + (-2 + a) b) )/(-2 + a)}, a] > q/8 &&
q > 0 && 0 < b < 1 &&
0 < s < 1 -
b && -Root[
q^3 - 2 b^2 q^3 + b^4 q^3 - 4 b q^3 s + 4 b^3 q^3 s -
2 q^3 s^2 + 6 b^2 q^3 s^2 + 4 b q^3 s^3 +
q^3 s^4 + (6 q^2 + 10 b^2 q^2 + 20 b q^2 s +
10 q^2 s^2) #1 + (12 q + b^2 q + 2 b q s + q s^2) #1^2 +
8 #1^3 &, 1] < (-1 + b + s)/(-1 + b), {b, 0, 1}, {s, 0,
1}], {{q, 1, "q"}, 0, 3}]

一个挺有意思的问题，虽然贴吧的编辑环境越来越不适合代码，还是试着在这答一下吧，希望楼不要被抽掉
首先注意到，求极值过程和q是没有关系的：由于q>0，可以作为常数提取到式子外，从而求极值的式子只包含参数b和s，这能很大程度上减少后续的计算量
然后处理一下a的范围。严格来讲，函数在有限开区间上的最值很多时候是不存在的，但实际计算中通常会把边界上的点也当作可以取到，所以这里直接把a的范围全部设置成闭区间，不会影响计算
接下来可以开始进行计算，为了简洁把函数和约束条件都写成函数的形式。
可以看到，无法满足约束条件的情况用时远大于正常情况，为了提高运算速度，需要在调用MaxValue
前排
除这种情况，也即计算出使a定义域不为空集的参数b、s需要满足的条件。这一步可以通过手算推导完成，这里介绍一种基于Reduce的简单方法如下，基本思路为将所有形如l<a<r的不等式变为l<r，再统一约化求解
可以看到，加入判断后避免了错误情况的耗时，但固定s绘制函数随b变化的曲线仍需耗时20秒左右，绘制RegionPlot仍不太现实，需要继续提速。
提速的重点放在MaxValue上，对于这种相对简单的问题，可以自己手写函数
手写myMaxValue2如上，可以看到其结果与myMaxValue一致，速度提升50倍左右。
注意到aStation所得结果为某6次多项式方程的6个根，使用NRoots替代Solve可以进一步提速3倍左右
开始画图，虽然还是要将近40秒，不过好歹画出来了
下层楼贴代码文本

Clear["`*"]
expr[b_, s_, a_] :=
1/2 a (1 - 2 (-1 + a) b^2 + (
a^2 b s)/(-1 +
a) - ((2 + (-2 + a) a) s^2)/(-1 + a)^2 + ((-1 + a) b)/(
b - a b + s));
cons[b_, s_, a_] :=
0 <= a <= 1 && 0 <= b <= -(1/(-2 + a)) &&
0 <= s <= ((-1 + a) (1 + (-2 + a) b))/(-2 + a);
bsLegalQ[b_, s_] =
Reduce[cons[b, s, a] /. LessEqual -> Less, a] /.
f_[l_, op_, a, op_, r_] :> f[l, op, r] // Reduce // FullSimplify;
(*expr对a的二阶导*)
ddexpr[b_, s_, a_] := Evaluate[FullSimplify@D[expr[b, s, a], a, a]];
(*驻点，即一阶导为0处*)
eqn[b_, s_] =
Evaluate[a /. Solve[D[expr[b, s, a], a] == 0, a]][[1, 1]]@a;
aStation[b_, s_] := List @@ (NRoots[eqn[b, s], a][[All, 2]]);
(*边界点*)
aBound[b_, s_] :=
Evaluate@DeleteDuplicates@
Flatten@Cases[Reduce[cons[b, s, a], a],
f_[l_, op_, a_, op_, r_] :> {l, r}, \[Infinity]];
myMaxValue[b_, s_] /; bsLegalQ[b, s] := Max[
expr[b, s,
(*选取满足约束条件且二阶导<0的驻点，与满足约束条件的边界点，计算函数值，所得列表最大值即为函数在区间内最大值*)
Join[Select[
aStation[b,
s], # \[Element] Reals && cons[b, s, #] &&
ddexpr[b, s, #] < 0 &],
Select[aBound[b, s], # \[Element] Reals && cons[b, s, #] &]]
]
]
\[ScriptCapitalR] = ImplicitRegion[bsLegalQ[b, s], {b, s}];
RegionPlot[
myMaxValue[b, s] > 1/8, {b, s} \[Element] \[ScriptCapitalR],
ImageSize -> 400] // AbsoluteTiming