differential equations - Simulating a combination of PDEs and ODEs

I am trying to simulate a combination of PDEs and ODEs, given below.

$$ \begin{matrix} -L\dfrac{\partial}{\partial t}I(t,z)&=&\dfrac{\partial}{\partial z}V(t,z)+RI(t,z)\\ C\dfrac{\partial}{\partial t}V(t,z)&=& -GV(t,z)-\dfrac{\partial}{\partial z}I(t,z) \end{matrix}\hspace{8mm}0

$$ \begin{matrix} -L_0\dfrac{d}{d t}I(t,0)&=&R_0I(t,0)+V(t,0)-E_0(t)\\ \dfrac{d}{dt}E_0(t)&=&-K_p\left(I(t,0)-I_d\right)-K_d \dfrac{\partial}{\partial t}I(t,0) \end{matrix}~\;\;\;z=0\\ $$

$$ \begin{matrix} -L_1\dfrac{d}{d t}I(t,1)&=& R_1I(t,1)-V(t,1)+v(t)\\ C_1\dfrac{d}{dt}v(t)&=& I(t,1)-\dfrac{v(t)}{R_C} \end{matrix}~\;\;\;\;z=1 $$

All constants can be taken as 1.

PS: I have tried the following:

l = 1; l0 = 1; l1 = 1; c = 1; c1 = 1;
r = 1; g = 1; r0 = 1; r1 = 1; rc = 1;
Ki = 1;
Kd = 2;
Ib0d = 5;
simtime = 10;
var = {Is[t, x], Vs[t, x], Ib0[t], Ib1[t], Vb1[t], Eb0[t]} // Flatten
varInit = var /. {t -> 0};

pde = Inverse[
    DiagonalMatrix[{-l, c}]].{D[Vs[t, x], x] + 
     r*Is[t, x], -g*Vs[t, x] - D[Is[t, x], x]};
MatrixForm[pde];
ode = Inverse[
    DiagonalMatrix[{-l0, -l1, c1}]].{r0*Ib0[t] + Vs[t, 0] - Eb0[t], 
    r1*Ib1[t] + Vb1[t] - Vs[t, 1], Ib1[t] - Vb1[t]/rc};
MatrixForm[ode];
uode = {-Ki*(Ib0[t] - Ib0d) - Kd*(r0*Ib0[t] + Vs[t, 0] - Eb0[t])/-l0};
MatrixForm[uode];

f = {pde, ode, uode} // Flatten;
DAE = {Thread[D[var, t] == f], Ib0[t] == Is[t, 0], Ib1[t] == Is[t, 1],
    Thread[ varInit == Flatten[RandomReal[{0, 1}, {1, 6}]]]} // 
  Flatten; MatrixForm[DAE]
uval = NDSolve[DAE, {x, 0, 1}, {t, 0, simtime}]

but it is not working, giving me an underdetermined system error.

Answer

This problem can be solved with the help of pdetoode.

We first eliminate Ib0 and Ib1 from the code because DAE system is troublesome compared to ODE system for NDSolve. It's better to avoid introducing them from the very beginning to make the code cleaner, but here I'll just modify your DAE for simplicity:

newdae = DeleteCases[
  DAE /. Rule @@@ (DAE[[7 ;; 8]] /. (a_ == b_) :> (Head@a == (Function[t, #] &@b))), 
  True | Is[__Integer] == _]

Then we extract those equations involving Is and Vs:

{eq@1, eq@2} = GatherBy[newdae, FreeQ[#, Is | Vs[__]] &]

Finally, discretize those involving Is and Vs and solve the resulting system:

domain = {0, 1}; points = 25; difforder = 2;
grid = Array[# &, points, domain];

ptoofunc = pdetoode[{Is, Vs}[t, x], t, grid, difforder];
delete = #[[2 ;; -2]] &;

(* Definition of pdetoode isn't included in this post,
   please find it in the link above. *)
eqlst = MapAt[delete, ptoofunc@eq@1, {1}]

uval = NDSolveValue[{eqlst, eq@2}, {Is /@ grid, Vs /@ grid, Head /@ var[[5 ;;]]}, {t, 0, 
   simtime}]

Is and Vs is discretized in $x$ direction in uval. If you prefer a continuous version, we can:

sol@Is = rebuild[uval[[1]], grid]

sol@Vs = rebuild[uval[[2]], grid]

Plot3D[sol[Vs][t, x], {t, 0, 10}, {x, 0, 1}]