87506e3d63d9b33636838ef6ee2bc9c3e3a0187e
[scilab.git] / scilab / modules / optimization / demos / icse / icsenb.f
1       subroutine icsenb(ind,nu,u,co,g,itv,rtv,dtv)
2 c     Copyright INRIA
3       external nbsec,nvcout,icsei
4 c     controle de la rentree de la navette en backward
5       call icse(ind,nu,u,co,g,itv,rtv,dtv,nbsec,nvcout,icsei)
6       end
7       subroutine icsenf(ind,nu,u,co,g,itv,rtv,dtv)
8       external nfsec,nvcout,icsei
9 c     controle de la rentree de la navette en forward
10       call icse(ind,nu,u,co,g,itv,rtv,dtv,nfsec,nvcout,icsei)
11       end
12
13       subroutine nbsec(indf,t,y,uc,uv,f,fy,fu,b,itu,dtu,
14      & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
15      & itmx,nex,nob,ntob,ntobi,nitu,ndtu)
16 c
17 c      Second membre de l'equation d'etat :
18 c       Parametres d'entree :
19 c        indf     : vaut 1,2,3 suivant qu'on veut calculer f,fy,fu
20 c        t        : instant courant
21 c        y(ny)    : etat a un instant donne
22 c        uc(nuc)  : controle independant du temps
23 c        uv(nuv)  : controle dependant du temps, a l'instant t
24 c        b(ny)    : terme constant dans le cas lineaire quadratique
25 c       Parametres de sortie :
26 c         indf    : >0 si  le calcul s'est  correctement effectue,<=0
27 c                   sinon
28 c        f(ny)    : second membre
29 c        fy(ny,ny): jacobien de f par rapport a y
30 c        fu(ny,nuc+nuv) : derivee de f par rapport au controle
31 c       Tableaux de travail reserves a l'utilisateur :
32 c        itu(nitu): tableau entier
33 c        dtu(ndtu): tableau double precision
34 c       (nitu et ndtu sont initialises par le common icsez).
35 c!
36       implicit double precision (a-h,o-z)
37       dimension y(ny),uc(*),uv(*),f(ny),fy(ny,ny),fu(ny,*),
38      &     b(ny),itu(*),dtu(*),iu(5)
39       call navetb(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
40      &  ny,nuc,nuv,nitu,ndtu)
41       end
42
43       subroutine nfsec(indf,t,y,uc,uv,f,fy,fu,b,itu,dtu,
44      & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
45      & itmx,nex,nob,ntob,ntobi,nitu,ndtu)
46 c
47 c      Second membre de l'equation d'etat :
48 c       Parametres d'entree :
49 c        indf     : vaut 1,2,3 suivant qu'on veut calculer f,fy,fu
50 c        t        : instant courant
51 c        y(ny)    : etat a un instant donne
52 c        uc(nuc)  : controle independant du temps
53 c        uv(nuv)  : controle dependant du temps, a l'instant t
54 c        b(ny)    : terme constant dans le cas lineaire quadratique
55 c       Parametres de sortie :
56 c         indf    : >0 si  le calcul s'est  correctement effectue,<=0
57 c                   sinon
58 c        f(ny)    : second membre
59 c        fy(ny,ny): jacobien de f par rapport a y
60 c        fu(ny,nuc+nuv) : derivee de f par rapport au controle
61 c       Tableaux de travail reserves a l'utilisateur :
62 c        itu(nitu): tableau entier
63 c        dtu(ndtu): tableau double precision
64 c       (nitu et ndtu sont initialises par le common icsez).
65 c!
66       implicit double precision (a-h,o-z)
67       dimension y(ny),uc(*),uv(*),f(ny),fy(ny,ny),fu(ny,*),
68      &     b(ny),itu(*),dtu(*),iu(5)
69       call navetf(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
70      &  ny,nuc,nuv,nitu,ndtu)
71       end
72
73       subroutine nvcout(indc,nu,tob,obs,cof,ytob,ob,u,
74      & c,cy,g,yob,d,itu,dtu,
75      & t0,tf,dti,dtf,ermx,iu,nuc,nuv,ilin,nti,ntf,ny,nea,
76      & itmx,nex,nob,ntob,ntobi,nitu,ndtu)
77 c
78 c
79 c     critere standard des moindres carres
80 c
81 c
82 c      Cout ponctuel :
83 c       Parametres d'entree :
84 c        indc     : 1 si on desire calculer c2,2 si on desire
85 c                   calculer c2y,c2u
86 c        tob      : instants de mesure
87 c        obs      : matrice d'observation
88 c        cof      : coefficients de ponderation du cout
89 c        ytob     : valeur de l'etat aux instants d'observation
90 c        ob       : mesures
91 c        u(nu)    : controle.Le controle variable est stocke a la
92 c                   suite du controle suite du constant.
93 c       Parametres de sortie :
94 c        indc     : comme pour icsec1
95 c        c2       : cout
96 c        c2y(ny,ntob) : derivee de c2 par rapport a y
97 c        g(nu)  : derivee de c2 par rapport a u
98 c!
99       implicit double precision (a-h,o-z)
100       dimension tob(ntob),obs(nob,ny),cof(nob,ntob),ytob(ny,ntob),
101      &ob(nex,ntob,nob),u(nu),cy(ny,ntob),g(nu),yob(nob,ntob),
102      &d(nob),itu(nitu),dtu(ndtu),iu(5)
103 c
104       call navetc(indc,nu,tob,obs,cof,ytob,ob,u,
105      &     c,cy,g,yob,d,itu,dtu)
106 c
107       end
108       subroutine navetb(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
109      &  ny,nuc,nuv,nitu,ndtu)
110 C        test de icse : navette
111 C        fb,inria
112 C!
113       implicit double precision (a-h,o-z)
114       dimension y(ny),uc(*),uv(*),f(ny),fy(ny,ny),fu(ny,nuc+nuv),itu(nit
115      &u),dtu(ndtu)
116       f(1)=-uc(1)*dtu(13)*(-dtu(8)*sin(y(2)*dtu(22))/(y(3)*dtu(23)+dtu(9
117      &))**2-0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu(2))*dtu
118      &(10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(21)
119       f(2)=-uc(1)*dtu(13)*(y(1)*dtu(21)*cos(y(2)*dtu(22))/(y(3)*dtu(23)+
120      &dtu(9))-dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21)*(y(3)*dtu(23)+dtu(9
121      &))**2)+0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*dtu(21)*exp(-
122      &y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(22)
123       f(3)=-uc(1)*y(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
124       f(4)=-uc(1)*y(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+
125      &dtu(9))*dtu(24))
126       if(indf.eq.2) then
127         fy(1,1)=dtu(1)*uc(1)*y(1)*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu(2))*
128      &  dtu(10)*dtu(13)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/dtu(7)
129         fy(1,2)=uc(1)*dtu(8)*dtu(13)*dtu(22)*cos(y(2)*dtu(22))/(dtu(21)*
130      &  (y(3)*dtu(23)+dtu(9))**2)
131         fy(1,3)=-uc(1)*dtu(13)*dtu(23)*(2*dtu(8)*sin(y(2)*dtu(22))/(y(3)
132      &  *dtu(23)+dtu(9))**3+0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dt
133      &  u(3)+dtu(2))*dtu(10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/(dtu(
134      &  7)*dtu(11)))/dtu(21)
135         fy(1,4)=0
136         fy(2,1)=-uc(1)*dtu(13)*dtu(21)*(cos(y(2)*dtu(22))/(y(3)*dtu(23)+
137      &  dtu(9))+dtu(8)*cos(y(2)*dtu(22))/(y(1)**2*dtu(21)**2*(y(3)*dtu(2
138      &  3)+dtu(9))**2)+0.5*dtu(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*exp(-y(3
139      &  )*dtu(23)/dtu(11))/dtu(7))/dtu(22)
140         fy(2,2)=-uc(1)*dtu(13)*(dtu(8)*sin(y(2)*dtu(22))/(y(1)*dtu(21)*(
141      &  y(3)*dtu(23)+dtu(9))**2)-y(1)*dtu(21)*sin(y(2)*dtu(22))/(y(3)*dt
142      &  u(23)+dtu(9)))
143         fy(2,3)=-uc(1)*dtu(13)*dtu(23)*(-y(1)*dtu(21)*cos(y(2)*dtu(22))/
144      &  (y(3)*dtu(23)+dtu(9))**2+2*dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21
145      &  )*(y(3)*dtu(23)+dtu(9))**3)-0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5)
146      &  )*dtu(10)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/(dtu(7)*dtu(11)))/d
147      &  tu(22)
148         fy(2,4)=0
149         fy(3,1)=-uc(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
150         fy(3,2)=-uc(1)*y(1)*dtu(13)*dtu(21)*dtu(22)*cos(y(2)*dtu(22))/dt
151      &  u(23)
152         fy(3,3)=0
153         fy(3,4)=0
154         fy(4,1)=-uc(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+
155      &  dtu(9))*dtu(24))
156         fy(4,2)=uc(1)*y(1)*dtu(13)*dtu(21)*dtu(22)*sin(y(2)*dtu(22))/((y
157      &  (3)*dtu(23)+dtu(9))*dtu(24))
158         fy(4,3)=uc(1)*y(1)*dtu(13)*dtu(21)*dtu(23)*cos(y(2)*dtu(22))/((y
159      &  (3)*dtu(23)+dtu(9))**2*dtu(24))
160         fy(4,4)=0
161       end if
162       if(indf.eq.3) then
163         fu(1,1)=-dtu(13)*(-dtu(8)*sin(y(2)*dtu(22))/(y(3)*dtu(23)+dtu(9)
164      &  )**2-0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu(2))*dt
165      &  u(10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(21)
166         fu(1,2)=0.5*dtu(1)*uc(1)*y(1)**2*(2*uv(1)*dtu(4)+dtu(3))*dtu(10)
167      &  *dtu(13)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/dtu(7)
168         fu(2,1)=-dtu(13)*(y(1)*dtu(21)*cos(y(2)*dtu(22))/(y(3)*dtu(23)+d
169      &  tu(9))-dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21)*(y(3)*dtu(23)+dtu(
170      &  9))**2)+0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*dtu(21)*ex
171      &  p(-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(22)
172         fu(2,2)=-0.5*dtu(1)*uc(1)*y(1)*dtu(6)*dtu(10)*dtu(13)*dtu(21)*ex
173      &  p(-y(3)*dtu(23)/dtu(11))/(dtu(7)*dtu(22))
174         fu(3,1)=-y(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
175         fu(3,2)=0
176         fu(4,1)=-y(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+d
177      &  tu(9))*dtu(24))
178         fu(4,2)=0
179       end if
180       end
181
182
183       subroutine navetf(indf,t,y,uc,uv,f,fy,fu,itu,dtu,
184      &  ny,nuc,nuv,nitu,ndtu)
185 C        test de icse : navette
186 C        fb,inria
187 C        [" "," "]
188 C
189       implicit double precision (a-h,o-z)
190       dimension y(ny),uc(*),uv(*),f(ny),fy(ny,ny),fu(ny,nuc+nuv),itu(nit
191      &u),dtu(ndtu)
192       f(1)=uc(1)*dtu(13)*(-dtu(8)*sin(y(2)*dtu(22))/(y(3)*dtu(23)+dtu(9)
193      &)**2-0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu(2))*dtu(
194      &10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(21)
195       f(2)=uc(1)*dtu(13)*(y(1)*dtu(21)*cos(y(2)*dtu(22))/(y(3)*dtu(23)+d
196      &tu(9))-dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21)*(y(3)*dtu(23)+dtu(9)
197      &)**2)+0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*dtu(21)*exp(-y
198      &(3)*dtu(23)/dtu(11))/dtu(7))/dtu(22)
199       f(3)=uc(1)*y(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
200       f(4)=uc(1)*y(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+d
201      &tu(9))*dtu(24))
202       if(indf.eq.2) then
203         fy(1,1)=-1.0*dtu(1)*uc(1)*y(1)*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu
204      &  (2))*dtu(10)*dtu(13)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/dtu(7)
205         fy(1,2)=-uc(1)*dtu(8)*dtu(13)*dtu(22)*cos(y(2)*dtu(22))/(dtu(21)
206      &  *(y(3)*dtu(23)+dtu(9))**2)
207         fy(1,3)=uc(1)*dtu(13)*dtu(23)*(2*dtu(8)*sin(y(2)*dtu(22))/(y(3)*
208      &  dtu(23)+dtu(9))**3+0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dtu
209      &  (3)+dtu(2))*dtu(10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/(dtu(7
210      &  )*dtu(11)))/dtu(21)
211         fy(1,4)=0
212         fy(2,1)=uc(1)*dtu(13)*dtu(21)*(cos(y(2)*dtu(22))/(y(3)*dtu(23)+d
213      &  tu(9))+dtu(8)*cos(y(2)*dtu(22))/(y(1)**2*dtu(21)**2*(y(3)*dtu(23
214      &  )+dtu(9))**2)+0.5*dtu(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*exp(-y(3)
215      &  *dtu(23)/dtu(11))/dtu(7))/dtu(22)
216         fy(2,2)=uc(1)*dtu(13)*(dtu(8)*sin(y(2)*dtu(22))/(y(1)*dtu(21)*(y
217      &  (3)*dtu(23)+dtu(9))**2)-y(1)*dtu(21)*sin(y(2)*dtu(22))/(y(3)*dtu
218      &  (23)+dtu(9)))
219         fy(2,3)=uc(1)*dtu(13)*dtu(23)*(-y(1)*dtu(21)*cos(y(2)*dtu(22))/(
220      &  y(3)*dtu(23)+dtu(9))**2+2*dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21)
221      &  *(y(3)*dtu(23)+dtu(9))**3)-0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5))
222      &  *dtu(10)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/(dtu(7)*dtu(11)))/dt
223      &  u(22)
224         fy(2,4)=0
225         fy(3,1)=uc(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
226         fy(3,2)=uc(1)*y(1)*dtu(13)*dtu(21)*dtu(22)*cos(y(2)*dtu(22))/dtu
227      &  (23)
228         fy(3,3)=0
229         fy(3,4)=0
230         fy(4,1)=uc(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+d
231      &  tu(9))*dtu(24))
232         fy(4,2)=-uc(1)*y(1)*dtu(13)*dtu(21)*dtu(22)*sin(y(2)*dtu(22))/((
233      &  y(3)*dtu(23)+dtu(9))*dtu(24))
234         fy(4,3)=-uc(1)*y(1)*dtu(13)*dtu(21)*dtu(23)*cos(y(2)*dtu(22))/((
235      &  y(3)*dtu(23)+dtu(9))**2*dtu(24))
236         fy(4,4)=0
237       end if
238       if(indf.eq.3) then
239         fu(1,1)=dtu(13)*(-dtu(8)*sin(y(2)*dtu(22))/(y(3)*dtu(23)+dtu(9))
240      &  **2-0.5*dtu(1)*y(1)**2*(uv(1)**2*dtu(4)+uv(1)*dtu(3)+dtu(2))*dtu
241      &  (10)*dtu(21)**2*exp(-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(21)
242         fu(1,2)=-0.5*dtu(1)*uc(1)*y(1)**2*(2*uv(1)*dtu(4)+dtu(3))*dtu(10
243      &  )*dtu(13)*dtu(21)*exp(-y(3)*dtu(23)/dtu(11))/dtu(7)
244         fu(2,1)=dtu(13)*(y(1)*dtu(21)*cos(y(2)*dtu(22))/(y(3)*dtu(23)+dt
245      &  u(9))-dtu(8)*cos(y(2)*dtu(22))/(y(1)*dtu(21)*(y(3)*dtu(23)+dtu(9
246      &  ))**2)+0.5*dtu(1)*y(1)*(uv(1)*dtu(6)+dtu(5))*dtu(10)*dtu(21)*exp
247      &  (-y(3)*dtu(23)/dtu(11))/dtu(7))/dtu(22)
248         fu(2,2)=0.5*dtu(1)*uc(1)*y(1)*dtu(6)*dtu(10)*dtu(13)*dtu(21)*exp
249      &  (-y(3)*dtu(23)/dtu(11))/(dtu(7)*dtu(22))
250         fu(3,1)=y(1)*dtu(13)*dtu(21)*sin(y(2)*dtu(22))/dtu(23)
251         fu(3,2)=0
252         fu(4,1)=y(1)*dtu(13)*dtu(21)*cos(y(2)*dtu(22))/((y(3)*dtu(23)+dt
253      &  u(9))*dtu(24))
254         fu(4,2)=0
255       end if
256       end
257
258       subroutine navetc(indc,nu,tob,obs,cof,ytob,ob,u,
259      & c,cy,g,yob,d,itu,dtu)
260 c
261 c     test probleme navette
262 c
263 c     sous programme appele par icse.f qui donne :
264 c     pour indc=1,le cout:c
265 c     pour indc=2,la matrice derivee du cout par rapport a l'etat
266 c                 calcule aux instants de mesure:cy(ny,ntob)
267 c     sorties :
268 c
269 c     c        double precision
270 c              cout
271 c     cy       double precision (ny,ntob)
272 c              derivee du cout par rapport a l'etat calcule aux
273 c              instants de mesure
274 c     g        double precision (nu)
275 c              le gradient est initialise a la derivee partielle du
276 c              cout par rapport au controle
277 c
278 c     variables internes:     yob,d
279 c
280       implicit double precision (a-h,o-z)
281       dimension ytob(4),ob(3),u(nu),cy(4),g(nu),dtu(*)
282 c
283       cpen=dtu(25)
284       cy(1)=ytob(1)-ob(1)
285       cy(2)=ytob(2)-ob(2)
286       cy(3)=ytob(3)-ob(3)
287       cy(4)=1.0d+0/cpen
288       c=ytob(4)/cpen + ( cy(1)**2 + cy(2)**2 + cy(3)**2 )/2.0d+0
289       do 10 k=1,nu
290 10    g(k)=0.0d+0
291       end