Formalization of one optimization problem or solution of inequalities

Question

I have polynomial:

$p=A_2t^2+A_1t+A_0$

$A_0=(x^2-y^2)+xz$

$A_1=x^2+y^2+z^2+\sin(x)$

$A_2=x^4+y^3+z^2$

$x,y,z$ - parameters, moreover $z$ - the value of which varies in the range $[0,1]$.

Polynomial $p$ has all real solutions if and only if:

$D=A_1^2-4A_0A_2>0$

There are known bounds:

$-1.25<A_0<2,-0.2<A_1<3.8,-1<A_2<3$ and $-1<x<1,-1<y<1$

Task: find the intervals for the parameters $x$ and $y$ at which the polynomial $p$ has only real solutions, i.e. the inequality $D>0$ (taking into account the existing bounds) for any $z$ from the specified range.

Problem: I tried to solve this problem numerically using the NMinimize and NMaximize, but to fill the intervals I need to restart the optimization many times, which I want to avoid.

How to solve the problem and get the required intervals for $x$ and $y$ out of the box? I.e. the answer should be something like:

x=[?;?] and y=[?;?]

My code:

Clear["Derivative"]
ClearAll["Global` "]

Remove[A, x, y, z]

Subscript[A, 0] = (x^2 - y^2) + x z;

Subscript[A, 1] = x^2 + y^2 + z^2 + Sin[x];

Subscript[A, 2] = x^4 + y^3 + z^2;

d = Subscript[A, 1]^2 - 4 Subscript[A, 0] Subscript[A, 2];

NMinimize[{x, 
   d > 0, -1.25 < Subscript[A, 0] < 2, -0.2 < Subscript[A, 1] < 
    3.8, -1 < Subscript[A, 2] < 3, -1 < x < 1, -1 < y < 1}, {x, y}, 
  Method -> {"RandomSearch", "SearchPoints" -> 5}];

NMaximize[{x, 
   d > 0, -1.25 < Subscript[A, 0] < 2, -0.2 < Subscript[A, 1] < 
    3.8, -1 < Subscript[A, 2] < 3, -1 < x < 1, -1 < y < 1}, {x, y}, 
  Method -> {"RandomSearch", "SearchPoints" -> 5}];

NMinimize[{y, 
   d > 0, -1.25 < Subscript[A, 0] < 2, -0.2 < Subscript[A, 1] < 
    3.8, -1 < Subscript[A, 2] < 3, -1 < x < 1, -1 < y < 1}, {x, y}, 
  Method -> {"RandomSearch", "SearchPoints" -> 5}];

NMaximize[{y, 
   d > 0, -1.25 < Subscript[A, 0] < 2, -0.2 < Subscript[A, 1] < 
    3.8, -1 < Subscript[A, 2] < 3, -1 < x < 1, -1 < y < 1}, {x, y}, 
  Method -> {"RandomSearch", "SearchPoints" -> 5}];

Answer 1

You can get a sense of the ranges just by plotting where the discriminant vanishes.

aa[0] = (x^2 - y^2) + x*z;
aa[1] = x^2 + y^2 + z^2 + Sin[z];
aa[2] = x^4 + y^3 + z^2;
poly = Array[aa, 3, 0].t^Range[0, 2];
disc = Discriminant[poly, t];
bounds = {-5/4 <= aa[0] <= 2,
  -1/5 <= aa[1] <= 19/5, -1 <= aa[2] <= 3,
  -1 <= x <= 1, -1 <= y <= 1, 0 <= z <= 1};

ContourPlot3D[disc == 0, {x, -1, 1},
  {y, -1, 1}, {z, 0, 1}]

I simple check at the origin shows that the larger central part is where the discriminant is positive.