1c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab 2c Copyright (C) INRIA 3c 4c Copyright (C) 2012 - 2016 - Scilab Enterprises 5c 6c This file is hereby licensed under the terms of the GNU GPL v2.0, 7c pursuant to article 5.3.4 of the CeCILL v.2.1. 8c This file was originally licensed under the terms of the CeCILL v2.1, 9c and continues to be available under such terms. 10c For more information, see the COPYING file which you should have received 11c along with this program. 12c 13C/MEMBR ADD NAME=ICSE2,SSI=0 14c 15 subroutine icse2(ind,nu,u,co,g,icsef,icsec2,icsei,y0,tob,obs, 16 &ob,ytot,f,b,fy,fu,ipv2,itob,cof,ytob,c2y,y0u,dmy,smy,oldmu, 17 &y,oldp,p,p0,gt,yob,d,itu,dtu, 18 &t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea, 19 &itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi) 20c 21c sous programme de icse.fortran:calcul du gradient par 22c integration du systeme adjoint 23c 24 implicit double precision (a-h,o-z) 25 dimension iu(5), u(nu),g(nu),y0(ny),tob(ntob),obs(nob,ny), 26 &ob(nex,ntob,nob),ytot(ny,nti+ntf),f(ny),b(ny),fy(ny,ny), 27 &fu(ny,nuc+nuv),ipv2(ny),itob(ntob),cof(nob,ntob),ytob(ny,ntob), 28 &c2y(ny,ntob),y0u(ny,nu),dmy(ny,ny),smy(ny,ny), 29 &oldmu(ny,nuc+nuv),y(ny),oldp(ny),p(ny),p0(ny), 30 >(nu),yob(nob,ntob),d(nob),itu(nitu),dtu(ndtu) 31 external icsef,icsec2,icsei 32c 33 character*6 nomf,nomc,nomi 34c 35c Initialisation 36c 37 call dset(nu,0.0d+0,g,1) 38 call dset(ny,0.0d+0,p,1) 39 kt=nti+ntf 40 ktob=ntob 41c lui et nui servent quand l'etat initial depend du controle 42 if (iu(2).gt.0) lui=min(nu,1+nuc) 43 if (iu(1).gt.0) lui=1 44 nui=iu(1)*nuc+iu(2)*nuv*(nti+ntf+1) 45c 46c 47c Calcul de itob,vecteur des indices des instants de mesure 48c a partir de tob,vecteur des instants de mesure 49c itob(j) est l'entier le plus proche de tob(j)/dt 50c 51 do 1 j=1,ntobi 521 itob(j)=int(((tob(j)-t0)/dti)+0.50d+0) 53 if (ntobi.lt.ntob) then 54 itob(ntobi+1)=nti+int(0.50d+0+(tob(ntobi+1)-t0-nti*dti)/dtf) 55 endif 56 if (ntobi+1.lt.ntob) then 57 do 2 j=ntobi+2,ntob 582 itob(j)=itob(ntobi+1)+int(0.50d+0+(tob(j)-tob(ntobi+1))/dtf) 59 endif 60c 61c Ecriture de ytob tableau des valeurs de l'etat 62c aux instants de mesure 63c 64 do 10 j=1,ntob 65 do 10 i=1,ny 6610 call dcopy(ny,ytot(1,itob(j)),1,ytob(1,j),1) 67c 68c Si ind=2,on calcule seulement le cout 69c Si ind=3,on calcule seulement le gradient 70c Si ind=4,on calcule le cout et le gradient 71c 72 if (ind.ne.3) then 73 indc=1 74 call icsec2(indc,nu,tob,obs,cof,ytob,ob,u, 75 & co,c2y,g,yob,d,itu,dtu, 76 & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea, 77 & itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi) 78 if (indc.le.0) then 79 ind=indc 80 return 81 endif 82 endif 83 if (ind.eq.2) return 84c 85c Calcul du gradient du cout en utilisant l'etat adjoint 86c discretise: 87c calcul de la derivee partielle c2y du cout par rapport 88c a l'etat 89c 90 indc=2 91 call icsec2(indc,nu,tob,obs,cof,ytob,ob,u,co,c2y,g,yob,d,itu, 92 &dtu, 93 & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea, 94 & itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi) 95 if (indc.le.0) then 96 ind=indc 97 return 98 endif 99c 100c +Evaluations successives de la contribution au gradient 101c a chaque pas de temps a l'aide de l'etat adjoint(non stocke) 102c 103 do 100 kt=nti+ntf,1,-1 104c 105c *Calcul de l'etat adjoint p_kt au pas kt 106c Calcul du second membre p 107c Initialisation: 108c nt=nti+ntf 109c si kt=nt,p est nul 110c si kt<nt,on a au depart p=p_kt+1,etat adjoint au pas kt+1; 111c on prend p=p0,ou p0=(smy)t*I*p avec smy=Id+(dt/2).dfy(t,y,u) 112c avec dt=t_(kt+1)-t_kt,y=y_kt,t=t_kt,I designant la matrice 113c diagonale d'ordre ny dont les nea premiers coefficients 114c valent 0 et les autres 1 et Id designant la matrice 115c identite d'ordre ny; 116c si le systeme est affine avec partie lineaire autonome 117c (ilin=2) smy est calculee seulement aux premiers pas de temps 118c pour dti et dtf,sinon (ilin=0 ou 1) smy est calculee a chaque 119c pas de temps 120c 121c stockage de l'etat adjoint et de dt/2 au pas kt+1 122c 123 call dcopy(ny,p,1,oldp,1) 124 luv=min(nu,1+nuc+kt*nuv) 125c 126c calcul de y=y_kt,dt=t_(kt+1)-t_kt,dt2=dt/2, 127c dt2new=(t_kt-t_(kt-1))/2 128c 129 call dcopy(ny,ytot(1,kt),1,y,1) 130c 131 if (kt.lt.nti) then 132 t=kt*dti+t0 133 dt=dti 134 else 135 t=nti*dti+(kt-nti)*dtf+t0 136 dt=dtf 137 endif 138 dt2=dt/2.0d+0 139 if (kt.ne.nti) then 140 dt2new=dt2 141 else 142 dt2new=dti/2.0d+0 143 endif 144c 145c Dans le cas ilin<=1, 146c calcul de fy=dfy(t,y,u) puis de smy=Id+(dt/2).dmy 147c lorsque kt<(nti+ntf) 148c Sinon (ilin>1),fy=dfy a ete calcule dans icse1 149c 150c 151 if (ilin.le.1) then 152 indf=2 153 call icsef(indf,t,y,u,u(luv),f,fy,fu,b,itu,dtu, 154 & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea, 155 & itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi) 156 if (indf.le.0) then 157 ind=indf 158 return 159 endif 160 endif 161c 162 if (kt.ne.nti+ntf) then 163 if (ilin.le.1.or.kt.eq.nti+ntf-1.or.kt.eq.nti-1) then 164 do 30 i=1,ny 165 do 30 j=1,ny 16630 smy(i,j)=dt2*fy(i,j) 167 do 35 i=1,ny 16835 smy(i,i)=1.0d+0+smy(i,i) 169 endif 170c 171c calcul de p0=(smy)t*I*p puis p=p0 172c 173 if (nea.gt.0) then 174 do 40 i=1,nea 17540 p(i)=0.0d+0 176 endif 177 call dmmul(p,1,smy,ny,p0,1,1,ny,ny) 178c 179 call dcopy(ny,p0,1,p,1) 180 endif 181c 182c Fin du calcul du second membre p 183c si kt=itob(ktob),on ajoute c2y(.,ktob) au second membre p 184c 185 if (ktob.gt.0) then 186 if (kt.eq.itob(ktob)) then 187 do 50 i=1,ny 18850 p(i)=p(i)+c2y(i,ktob) 189 ktob=ktob-1 190 endif 191 endif 192c 193c Calcul et factorisation de la matrice dmy du systeme 194c de l'etat adjoint 195c dmy=I-dt2new.dfy(t,y,u) avec dt2new=(t_kt-t_(kt-1))/2, 196c y=y_kt,t=t_kt,Idesignant la matrice diagonale d'ordre 197c ny dont les nea premiers coefficients valent 0 et les 198c autres 1; 199c si le systeme est affine avec partie lineaire autonome 200c (ilin=2) dmy est calculee et factorisee aux premiers 201c pas de temps pour dti et dtf,sinon (ilin=0 ou 1) dmy est 202c calculee et factorisee a chaque pas de temps 203c 204 if (ilin.le.1.or.kt.eq.nti+ntf.or.kt.eq.nti) then 205 do 60 i=1,ny 206 do 60 j=1,ny 20760 dmy(i,j)=-dt2new*fy(i,j) 208 do 65 i=1+nea,ny 20965 dmy(i,i)=1.0d+0+dmy(i,i) 210 call dgefa(dmy,ny,ny,ipv2,info) 211 endif 212c 213c Resolution de (dmy)t*X=p,la solution s'appelant p 214c p est alors p_kt,etat adjoint au pas kt 215c 216 call dgesl(dmy,ny,ny,ipv2,p,1) 217c 218c *Calcul du gradient g au pas kt+1 219c calcul de la contribution gt au gradient au pas kt+1: 220c gt=(dt/2).(I*dfu(t_kt,y_kt,u)+dfu(t_kt+1,y_kt+1,u))t*p_kt+1 221c avec dt=t_(kt+1)-t_kt,I designant la matrice diagonale 222c d'ordre ny dont les nea premiers coefficients valent 0 223c et les autres 1 224c 225 if (nuv.gt.0.or.iu(3).eq.1) then 226 indf=3 227 call icsef(indf,t,y,u,u(luv),f,fy,fu,b,itu,dtu, 228 & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea, 229 & itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi) 230 if (indf.le.0) then 231 ind=indf 232 return 233 endif 234 if (kt.lt.nti+ntf) then 235 call dmmul(oldp,1,oldmu,ny,gt,1,1,ny,nuc+nuv) 236 call dscal(nuc+nuv,dt2,gt,1) 237c le gradient g est la somme des contributions 238 if(iu(3).gt.0) then 239 call dadd(nuc,gt,1,g,1) 240 endif 241 if(nuv.gt.0) then 242 luv=min(nu,1+nuc+(kt+1)*nuv) 243 call dadd(nuv,gt(nuc+1),1,g(luv),1) 244 endif 245 if (nea.gt.0) then 246 do 70 i=1,nea 24770 oldp(i)=0.0d+0 248 endif 249 call dmmul(oldp,1,fu,ny,gt,1,1,ny,nuc+nuv) 250 call dscal(nuc+nuv,dt2,gt,1) 251c le gradient g est la somme des contributions 252 if(iu(3).gt.0) then 253 call dadd(nuc,gt,1,g,1) 254 endif 255 if(nuv.gt.0) then 256 luv=min(nu,1+nuc+kt*nuv) 257 call dadd(nuv,gt(nuc+1),1,g(luv),1) 258 endif 259 endif 260c 261c stockage de dfu(t_kt,y_kt,u) dans oldmu 262c 263 call dcopy(ny*(nuc+nuv),fu,1,oldmu,1) 264c 265c *On passe au pas de temps suivant:kt-1,sauf si kt=1,auquel cas 266c on calcule la contribution gt au gradient au pas kt=1 et 267c on l'ajoute a g;on a: 268c gt=(dt/2).(I*dfu(t0,y0,u)+dfu(t_1,y_1,u))t*p_1,avec 269c dt=t_1-t0=dti car nti n'est jamais nul par convention 270c 271 if (kt.eq.1) then 272 t=t0 273 dt2=dti/2.0d+0 274 call dcopy(ny,y0,1,y,1) 275 indf=3 276 call icsef(indf,t,y,u,u(luv),f,fy,fu,b,itu,dtu, 277 & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea, 278 & itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi) 279 if (indf.le.0) then 280 ind=indf 281 return 282 endif 283 call dmmul(p,1,oldmu,ny,gt,1,1,ny,nuc+nuv) 284 call dscal(nuc+nuv,dt2,gt,1) 285c le gradient g est la somme des contributions 286 if(iu(3).gt.0) then 287 call dadd(nuc,gt,1,g,1) 288 endif 289 if(nuv.gt.0) then 290 luv=min(nu,1+nuc+nuv) 291 call dadd(nuv,gt(nuc+1),1,g(luv),1) 292 endif 293 if (nea.gt.0) then 294 do 90 i=1,nea 29590 p(i)=0.0d+0 296 endif 297 call dmmul(p,1,fu,ny,gt,1,1,ny,nuc+nuv) 298 call dscal(nuc+nuv,dt2,gt,1) 299c le gradient g est la somme des contributions 300 if(iu(3).gt.0) then 301 call dadd(nuc,gt,1,g,1) 302 endif 303 if(nuv.gt.0) then 304 luv=min(nu,1+nuc) 305 call dadd(nuv,gt(nuc+1),1,g(luv),1) 306 endif 307 endif 308 endif 309100 continue 310c 311c gradient par rapport au controle initial 312c 313 if(max(iu(1),iu(2)) .gt.0) then 314c 315c calcul de l'etat adjoint initial p0 316c 317 indf=2 318 call icsef(indf,t,y,u,u(luv),f,fy,fu,b,itu,dtu, 319 & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea, 320 & itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi) 321 if (indf.eq.0) then 322 ind=indf 323 return 324 endif 325 do 120 i=1,ny 326 do 120 j=1,ny 327120 smy(i,j)=dt2*fy(i,j) 328 do 125 i=1,ny 329125 smy(i,i)=1.0d+0+smy(i,i) 330 if (nea.gt.0) then 331 do 130 i=1,nea 332130 p(i)=0.0d+0 333 endif 334 call dmmul(p,1,smy,ny,p0,1,1,ny,ny) 335c incrementation du gradient 336 indi=2 337 call icsei(indi,nui,u(lui),y0,y0u,itu,dtu, 338 & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea, 339 & itmx,nex,nob,ntob,ntobi,nitu,ndtu,nomf,nomc,nomi) 340 if (indi.le.0) then 341 ind=indi 342 return 343 endif 344 call dmmul(p0,1,y0u,ny,gt,1,1,nui,nui) 345 do 150 i=1,nui 346150 g(lui+i-1)=g(lui+i-1)+gt(i) 347c 348 endif 349 end 350