ODE coupled to timedependent polynomial

First we have to decide which initial conditions to use for the DAE system.

soly=Values@Solve[ (-82.6333 + 0.3/a[t]^3) y[t] + 50 y[t]^2 + 49.3 y[t]^3 - (50 y[t]^4)/3 == 0, y[t]] /. a[t] -> 1 // Chop
(*{{0}, {-1.45186}, {3.41291}, {0.996957}}*)

OP looks for a branch 0<y<1 (see comment)

y1=soly[[-1]][[1]] (*0.996957*)

Now we solve the DAE-system with intial conditions a[1]==1, y[1]==y1

sola = NDSolveValue[{-0.7 - 0.3/a[t]^3 -50/3 (-1 + y[t]) (2 + (2 - y[t]) y[t]) + a'[t]^2/a[t]^2 ==0
, (-82.6333 + 0.3/a[t]^3) y[t] + 50 y[t]^2 +49.3 y[t]^3 - (50 y[t]^4)/3 == 0 , a[1] == 1, y[1] == y1},a, {t, 1, 2}]
Plot[sola[t] , {t, 1, 2}]

Note: I am not very experienced with NDSolve and differential algebraic equations, so the following solution is probably not the nicest one.

---

Normally, NDSolve is able to handle differential algebraic equations automatically. However, such equations usually have more than one solution, and we need to do some manual work to select the one we need. While the below works, I am also interested in the best way to do this.

Taking your equations,

deq = -0.7 - 0.3/a[t]^3 - 50/3 (-1 + y[t]) (2 + (2 - y[t]) y[t]) + 
   a'[t]^2/a[t]^2 == 0

eq = (-82.6333 + 0.3/a[t]^3) y[t] + 50 y[t]^2 + 
   49.3 y[t]^3 - (50 y[t]^4)/3 == 0

we can try a direct solution:

sol = NDSolve[
  {deq, eq, a[0] == 1},
  {a, y},
  {t, 0, 10}
  ]

But there is a warning that the computed initial conditions do not match the initial conditions we provided manually:

Let's verify them:

In[4]:= {a[0], y[0]} /. sol
Out[4]= {{2.7121*10^6, 1.}}

Indeed, a[0] is far from 1. I assume (but I don't know) that this might be due to Mathematica choosing the "wrong" y[0] solution to eq.

Let us do some manual investigation then. eq is a fourth-order polynomial, so for a given a[0] we have four possible y[0] values. But not all real y[0] values are acceptable, as some of them lead to a negative a'[0] when solving deq.

This plot shows the range of acceptable y[0] values:

Plot[
 ad /. Solve[deq /. t -> 0 /. a[0] -> 1 /. a'[0]^2 -> ad, ad] // 
  Evaluate,
 {y[0], -3, 4}
 ]

For a[0] == 1, we get these possible y[0] (all real):

ysol = Solve[eq /. t -> 0 /. a[0] -> 1]
(* {{y[0] -> -1.45186}, {y[0] -> 0.}, {y[0] -> 0.996957}, {y[0] -> 3.41291}} *)

But only the first and third one give valid a'[0] values:

NSolve /@ (deq /. t -> 0 /. a[0] -> 1 /. ysol)
(* {{{Derivative[1][a][0] -> -11.1386}, {Derivative[1][a][0] -> 
    11.1386}}, {{Derivative[1][a][0] -> 
    0. - 5.68624 I}, {Derivative[1][a][0] -> 
    0. + 5.68624 I}}, {{Derivative[1][a][
     0] -> -0.920791}, {Derivative[1][a][0] -> 
    0.920791}}, {{Derivative[1][a][0] -> 
    0. - 10.6062 I}, {Derivative[1][a][0] -> 0. + 10.6062 I}}} *)

We can compute a solution e.g. using the third one:

NDSolve[
 {deq, eq, a[0] == 1, ysol[[3]] /. Rule -> Equal},
 {a, y},
 {t, 0, 10}
 ]

(EDIT: I originally missed OP's comment that y[0] should be in [0,1]. It seems that it is indeed the 3rd solution that satisfies this.)

Note that this is still not the only solution for the given y[0], as a'[0] can be positive or negative. You can compute a desired a'[0] and add that to the initial conditions as well.

---

Now my personal question to other people who might write an answer:

Is there a simpler way to choose our desired solution? Is there a way to avoid manually computing the initial conditions and plugging them in? For example, can we somehow just specify that we want a'[0] >= 0 instead of computing an explicit value for a'[0] by hand and plugging that in? Can we do the same asking that 0 < y[0] < 1?