← Back to branch summary

~verifypn-cpn/verifypn/colored

« back to all changes in this revision

Viewing changes to PetriEngine/Structures/StateConstraints.cpp

  • Committer: Jonas Finnemann Jensen
  • Date: 2011-09-15 13:30:00 UTC
  • Revision ID: jopsen@gmail.com-20110915133000-wnywm1odf82emiuw
Import of sources from github

Show diffs side-by-side

added added

removed removed

Lines of Context:
PetriEngine/Structures/StateConstraints.cpp
 
1
/* PeTe - Petri Engine exTremE
 
2
 * Copyright (C) 2011  Jonas Finnemann Jensen <jopsen@gmail.com>,
 
3
 *                     Thomas Søndersø Nielsen <primogens@gmail.com>,
 
4
 *                     Lars Kærlund Østergaard <larsko@gmail.com>,
 
5
 * 
 
6
 * This program is free software: you can redistribute it and/or modify
 
7
 * it under the terms of the GNU General Public License as published by
 
8
 * the Free Software Foundation, either version 3 of the License, or
 
9
 * (at your option) any later version.
 
10
 * 
 
11
 * This program is distributed in the hope that it will be useful,
 
12
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
 * GNU General Public License for more details.
 
15
 * 
 
16
 * You should have received a copy of the GNU General Public License
 
17
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
 
18
 */
 
19
#include "StateConstraints.h"
 
20
 
 
21
#include <lpsolve/lp_lib.h>
 
22
#include <assert.h>
 
23
#include <stdio.h>
 
24
 
 
25
/** Don't give lp_solve more than 30s to handle the problem */
 
26
#define OVER_APPROX_TIMEOUT                             30
 
27
/** Number of M-traps to account for each place */
 
28
#define OVER_APPROX_TRAP_FACTOR                 2
 
29
 
 
30
namespace PetriEngine{
 
31
namespace Structures{
 
32
 
 
33
/** Attempts to merge two StateConstraints, return NULL if can't be merged */
 
34
StateConstraints* StateConstraints::merge(const StateConstraints* c) const{
 
35
        StateConstraints* nc = new StateConstraints(nPlaces, nVars);
 
36
        for(size_t p = 0; p < nPlaces; p++){
 
37
                if(!nc->setPlaceMax(p, pcs[p].max) ||
 
38
                   !nc->setPlaceMax(p, c->pcs[p].max) ||
 
39
                   !nc->setPlaceMin(p, pcs[p].min) ||
 
40
                   !nc->setPlaceMin(p, c->pcs[p].min)){
 
41
                        delete nc;
 
42
                        return NULL;
 
43
                }
 
44
        }
 
45
        for(size_t v = 0; v < nVars; v++){
 
46
                if(!nc->setVarMax(v, vcs[v].max) ||
 
47
                   !nc->setVarMax(v, c->vcs[v].max) ||
 
48
                   !nc->setVarMin(v, vcs[v].min) ||
 
49
                   !nc->setVarMin(v, c->vcs[v].min)){
 
50
                        delete nc;
 
51
                        return NULL;
 
52
                }
 
53
        }
 
54
        return nc;
 
55
}
 
56
 
 
57
/** True, if every marking satisfied by subset is also satisfied by this */
 
58
bool StateConstraints::isSuperSet(const StateConstraints* subset) const{
 
59
        for(size_t p = 0; p < nPlaces; p++){
 
60
                if(pcs[p].max < subset->pcs[p].max)
 
61
                        return false;
 
62
                if(pcs[p].min > subset->pcs[p].min)
 
63
                        return false;
 
64
        }
 
65
        for(size_t v = 0; v < nVars; v++){
 
66
                if(vcs[v].max < subset->vcs[v].max)
 
67
                        return false;
 
68
                if(vcs[v].min > subset->vcs[v].min)
 
69
                        return false;
 
70
        }
 
71
        return true;
 
72
}
 
73
 
 
74
void StateConstraints::reset(){
 
75
        for(size_t p = 0; p < nPlaces; p++){
 
76
                pcs[p].max = CONSTRAINT_INFTY;
 
77
                pcs[p].min = 0;
 
78
        }
 
79
        for(size_t v = 0; v < nVars; v++){
 
80
                vcs[v].max = CONSTRAINT_INFTY;
 
81
                vcs[v].min = 0;
 
82
        }
 
83
}
 
84
 
 
85
/** Merge the two sets of StateConstraints such that one from A and one from B is always satisfied, when one in the return value is
 
86
 * @remarks This will take ownership of the provided StateConstraints, and delete them or own them
 
87
 */
 
88
std::vector<StateConstraints*> StateConstraints::mergeAnd(const std::vector<StateConstraints*>& A, const std::vector<StateConstraints*>& B){
 
89
        std::vector<StateConstraints*> retval;
 
90
        typedef std::vector<StateConstraints*>::const_iterator iter;
 
91
        for(iter a = A.begin(); a != A.end(); a++){
 
92
                for(iter b = B.begin(); b != B.end(); b++){
 
93
                        StateConstraints* nc = (*a)->merge(*b);
 
94
                        if(nc)
 
95
                                retval.push_back(nc);
 
96
                }
 
97
                delete *a;
 
98
        }
 
99
        for(iter b = B.begin(); b != B.end(); b++)
 
100
                delete *b;
 
101
        return retval;
 
102
}
 
103
 
 
104
/** Merge the two sets of StateConstraints such that either on from A or B is always satisfied, when one in the return value is
 
105
 * @remarks This will take ownership of the provided StateConstraints, and delete them or own them
 
106
 */
 
107
std::vector<StateConstraints*> StateConstraints::mergeOr(const std::vector<StateConstraints*>& A, const std::vector<StateConstraints*>& B){
 
108
        std::vector<StateConstraints*> retval, excluded;
 
109
        typedef std::vector<StateConstraints*>::const_iterator iter;
 
110
        for(iter a = A.begin(); a != A.end(); a++){
 
111
                bool exclude = false;
 
112
                for(iter b = B.begin(); b != B.end(); b++){
 
113
                        exclude |= (*b)->isSuperSet(*a) && !(*b)->isEqual(*a);
 
114
                        if(exclude) break;
 
115
                }
 
116
                if(!exclude)
 
117
                        retval.push_back(*a);
 
118
                else
 
119
                        excluded.push_back(*a);
 
120
        }
 
121
        for(iter b = B.begin(); b != B.end(); b++){
 
122
                bool exclude = false;
 
123
                for(iter a = A.begin(); a != A.end(); a++){
 
124
                        exclude |= (*a)->isSuperSet(*b);
 
125
                        if(exclude) break;
 
126
                }
 
127
                if(!exclude)
 
128
                        retval.push_back(*b);
 
129
                else
 
130
                        excluded.push_back(*b);
 
131
        }
 
132
        for(iter e = excluded.begin(); e != excluded.end(); e++)
 
133
                delete *e;
 
134
        return retval;
 
135
}
 
136
 
 
137
StateConstraints* StateConstraints::requirePlaceMin(const PetriNet& net, int place, MarkVal value){
 
138
        StateConstraints* c = new StateConstraints(net);
 
139
        c->setPlaceMin(place, value);
 
140
        return c;
 
141
}
 
142
 
 
143
StateConstraints* StateConstraints::requirePlaceMax(const PetriNet& net, int place, MarkVal value){
 
144
        StateConstraints* c = new StateConstraints(net);
 
145
        c->setPlaceMax(place, value);
 
146
        return c;
 
147
}
 
148
 
 
149
StateConstraints* StateConstraints::requirePlaceEqual(const PetriNet& net, int place, MarkVal value){
 
150
        StateConstraints* c = new StateConstraints(net);
 
151
        c->setPlaceMax(place, value);
 
152
        c->setPlaceMin(place, value);
 
153
        return c;
 
154
}
 
155
 
 
156
 
 
157
StateConstraints* StateConstraints::requireVarMin(const PetriNet& net, int var, MarkVal value){
 
158
        StateConstraints* c = new StateConstraints(net);
 
159
        c->setVarMin(var, value);
 
160
        return c;
 
161
}
 
162
 
 
163
StateConstraints* StateConstraints::requireVarMax(const PetriNet& net, int var, MarkVal value){
 
164
        StateConstraints* c = new StateConstraints(net);
 
165
        c->setVarMax(var, value);
 
166
        return c;
 
167
}
 
168
 
 
169
StateConstraints* StateConstraints::requireVarEqual(const PetriNet& net, int var, MarkVal value){
 
170
        StateConstraints* c = new StateConstraints(net);
 
171
        c->setVarMax(var, value);
 
172
        c->setVarMin(var, value);
 
173
        return c;
 
174
}
 
175
 
 
176
/** Attempts to solve using lp_solve, returns True if the net cannot satisfy these constraints! */
 
177
bool StateConstraints::isImpossible(const PetriNet& net,
 
178
                                                                        const MarkVal* m0,
 
179
                                                                        const VarVal*) const{
 
180
        assert(nPlaces == net.numberOfPlaces());
 
181
        assert(nVars == net.numberOfVariables());
 
182
 
 
183
        /*printf("-------\n");
 
184
        for(int p = 0; p < nPlaces; p++)
 
185
                printf("%s: [%i, %i]\n", net.placeNames()[p].c_str(), pcs[p].min, pcs[p].max);
 
186
        printf("-------\n");*/
 
187
 
 
188
        // Create linary problem
 
189
        lprec* lp;
 
190
        lp = make_lp(0, net.numberOfTransitions());     // One variable for each entry in the firing vector
 
191
        assert(lp);
 
192
        if(!lp) return false;
 
193
 
 
194
        // Set verbosity
 
195
        set_verbose(lp, IMPORTANT);
 
196
 
 
197
        // Set transition names (not strictly needed)
 
198
        for(size_t i = 0; i < net.numberOfTransitions(); i++)
 
199
                set_col_name(lp, i+1, const_cast<char*>(net.transitionNames()[i].c_str()));
 
200
 
 
201
        // Start adding rows
 
202
        set_add_rowmode(lp, TRUE);
 
203
 
 
204
        REAL row[net.numberOfTransitions() + 1];
 
205
        for(size_t p = 0; p < nPlaces; p++){
 
206
                // Set row zero
 
207
                memset(row, 0, sizeof(REAL) * net.numberOfTransitions() + 1);
 
208
                for(size_t t = 0; t < net.numberOfTransitions(); t++){
 
209
                        int d = net.outArc(t, p) - net.inArc(p, t);
 
210
                        row[1+t] = d;
 
211
                }
 
212
 
 
213
                if(pcs[p].min == pcs[p].max &&
 
214
                   pcs[p].max != CONSTRAINT_INFTY){
 
215
                        double target = pcs[p].min - m0[p];
 
216
                        add_constraint(lp, row, EQ,  target);
 
217
                }else{
 
218
                        // There's always a min, even zero is interesting
 
219
                        double target = pcs[p].min - m0[p];
 
220
                        add_constraint(lp, row, GE,  target);
 
221
                        if(pcs[p].max != CONSTRAINT_INFTY){
 
222
                                double target = pcs[p].max - m0[p];
 
223
                                add_constraint(lp, row, LE,  target);
 
224
                        }
 
225
                }
 
226
        }
 
227
 
 
228
        // Finished adding rows
 
229
        set_add_rowmode(lp, FALSE);
 
230
 
 
231
        // Create objective
 
232
        memset(row, 0, sizeof(REAL) * net.numberOfTransitions() + 1);
 
233
        for(size_t t = 0; t < net.numberOfTransitions(); t++)
 
234
                row[1+t] = 1;   // The sum the components in the firing vector
 
235
 
 
236
        // Set objective
 
237
        set_obj_fn(lp, row);
 
238
 
 
239
        // Minimize the objective
 
240
        set_minim(lp);
 
241
 
 
242
        // Set variables as integer variables
 
243
        for(size_t i = 0; i < net.numberOfTransitions(); i++)
 
244
                set_int(lp, 1+i, TRUE);
 
245
 
 
246
        // Write it to stdout
 
247
        /*printf("--------------\n");
 
248
        write_LP(lp, stdout);
 
249
        printf("--------------\n");*/
 
250
 
 
251
        //Set timeout, and handle a timeout
 
252
        set_timeout(lp, OVER_APPROX_TIMEOUT);
 
253
 
 
254
        // Attempt to solve the problem
 
255
        int result = solve(lp);
 
256
 
 
257
        // Limit on traps to test
 
258
        size_t traplimit = nPlaces * OVER_APPROX_TRAP_FACTOR;
 
259
        // Try to add a minimal trap constraint
 
260
        while((result == OPTIMAL || result == SUBOPTIMAL) && traplimit-- < 0){
 
261
                memset(row, 0, sizeof(REAL) * net.numberOfTransitions() + 1);
 
262
                // Get the firing vector
 
263
                get_variables(lp, row);
 
264
                // Compute the resulting marking
 
265
                MarkVal rMark[net.numberOfPlaces()];
 
266
                for(size_t p = 0; p < nPlaces; p++){
 
267
                        rMark[p] = m0[p];
 
268
                        for(size_t t = 0; t < net.numberOfTransitions(); t++)
 
269
                                rMark[p] += (net.outArc(t, p) - net.inArc(p, t)) * (int)row[t];
 
270
                }
 
271
 
 
272
                // Find an M-trap
 
273
                BitField trap(minimalTrap(net, m0, rMark));
 
274
 
 
275
                /*printf("Initial marking:\n");
 
276
                for(size_t p = 0; p < nPlaces; p++)
 
277
                        printf("\t %s = %i\n", net.placeNames()[p].c_str(), m0[p]);
 
278
                // Debug information
 
279
                printf("Firing vector:\n");
 
280
                for(size_t t = 0; t < net.numberOfTransitions(); t++)
 
281
                        printf("\t %s = %i\n", net.transitionNames()[t].c_str(), (int)row[t]);
 
282
                printf("Result marking:\n");
 
283
                for(size_t p = 0; p < nPlaces; p++)
 
284
                        printf("\t %s = %i\n", net.placeNames()[p].c_str(), rMark[p]);
 
285
                printf("Got trap: ");
 
286
                for(size_t p = 0; p < nPlaces; p++){
 
287
                        if(trap.test(p))
 
288
                                printf("%s, ", net.placeNames()[p].c_str());
 
289
                }
 
290
                printf("\n");*/
 
291
 
 
292
                //Break if there's no trap
 
293
                if(trap.none()) break;
 
294
 
 
295
                // Compute the new equation
 
296
                for(size_t t = 0; t < net.numberOfTransitions(); t++){
 
297
                        row[1+t] = 0;
 
298
                        for(size_t p = 0; p < nPlaces; p++)
 
299
                                if(trap.test(p))
 
300
                                        row[1+t] += net.outArc(t, p) - net.inArc(p, t);
 
301
                }
 
302
 
 
303
                // Add a new row with target as greater than equal to 1
 
304
                set_add_rowmode(lp, TRUE);
 
305
                add_constraint(lp, row, GE,  1);
 
306
                set_add_rowmode(lp, FALSE);
 
307
 
 
308
                // Attempt to solve the again
 
309
                result = solve(lp);
 
310
        }
 
311
 
 
312
        // Delete the linear problem
 
313
        delete_lp(lp);
 
314
        lp = NULL;
 
315
 
 
316
        // Return true, if it was infeasible
 
317
        return result == INFEASIBLE;
 
318
}
 
319
 
 
320
BitField StateConstraints::maxTrap(const PetriNet& net,
 
321
                                                                   BitField places,
 
322
                                                                   const MarkVal* resultMarking) const{
 
323
        BitField next(places.size());
 
324
        BitField prev(places);
 
325
        do{
 
326
                prev = places;
 
327
                for(size_t i = 0; i < places.size(); i++)
 
328
                        next.set(i, isInMaxTrap(net, i, places, resultMarking));
 
329
                places = next;
 
330
        }while(prev != next);
 
331
        return places;
 
332
}
 
333
 
 
334
bool StateConstraints::isInMaxTrap(const PetriNet& net,
 
335
                                                                   size_t place,
 
336
                                                                   const BitField& places,
 
337
                                                                   const MarkVal* resultMarking) const{
 
338
        if(!places.test(place))
 
339
                return false;
 
340
        /*
 
341
                0 if M(p_i) = 1
 
342
                0 if there is (p_i , t) ∈ F such that x_j = 0
 
343
                        for every p_j ∈ t•
 
344
                1 otherwise
 
345
        */
 
346
        if(resultMarking[place] > 0)
 
347
                return false;
 
348
 
 
349
        for(unsigned int t = 0; t < net.numberOfTransitions(); t++){
 
350
                if(net.inArc(place, t) > 0){
 
351
                        bool exclude = true;
 
352
                        for(unsigned int j = 0; j < net.numberOfPlaces(); j++){
 
353
                                if(net.outArc(t, j) > 0){
 
354
                                        exclude &= !places.test(j);
 
355
                                }
 
356
                        }
 
357
                        if(exclude)
 
358
                                return false;
 
359
                }
 
360
        }
 
361
        return true;
 
362
}
 
363
 
 
364
BitField StateConstraints::minimalTrap(const PetriNet& net,
 
365
                                                                           const MarkVal* marking,
 
366
                                                                           const MarkVal* resultMarking) const{
 
367
        const size_t nPlaces = net.numberOfPlaces();
 
368
        BitField trap(nPlaces);
 
369
        trap.set();
 
370
        trap = maxTrap(net, trap, resultMarking);
 
371
        if(!isMarked(trap, marking))
 
372
                return BitField(nPlaces).clear();
 
373
 
 
374
        //Get the exclusion candidates
 
375
        BitField EC(trap);
 
376
        BitField tmp(nPlaces);
 
377
        while(EC.any()){
 
378
                int exclude = EC.first();
 
379
                tmp = trap;
 
380
                tmp.clear(exclude);
 
381
                EC.clear(exclude);
 
382
                tmp = maxTrap(net, tmp, resultMarking);
 
383
                if(isMarked(tmp, marking)){
 
384
                        trap = tmp;
 
385
                        EC = tmp;
 
386
                }
 
387
        }
 
388
        return trap;
 
389
}
 
390
 
 
391
bool StateConstraints::isMarked(const BitField& places, const MarkVal* marking) const{
 
392
        for(size_t p = 0; p < places.size(); p++)
 
393
                if(places.test(p) && marking[p] > 0)
 
394
                        return true;
 
395
        return false;
 
396
}
 
397
 
 
398
bool StateConstraints::isSpecific() const{
 
399
        for(size_t p = 0; p < nPlaces; p++){
 
400
                if(pcs[p].max != pcs[p].min)
 
401
                        return false;
 
402
        }
 
403
        for(size_t v = 0; v < nVars; v++){
 
404
                if(vcs[v].max != vcs[v].min)
 
405
                        return false;
 
406
        }
 
407
        return true;
 
408
}
 
409
 
 
410
int StateConstraints::fireVectorSize(const PetriNet& net,
 
411
                                                                         const MarkVal* m0,
 
412
                                                                         const VarVal*) const{
 
413
        assert(nPlaces == net.numberOfPlaces());
 
414
        assert(nVars == net.numberOfVariables());
 
415
 
 
416
        // Create linary problem
 
417
        lprec* lp;
 
418
        lp = make_lp(0, net.numberOfTransitions());     // One variable for each entry in the firing vector
 
419
        assert(lp);
 
420
        if(!lp) return false;
 
421
 
 
422
        // Set verbosity
 
423
        set_verbose(lp, IMPORTANT);
 
424
 
 
425
        // Set transition names (not strictly needed)
 
426
        for(size_t i = 0; i < net.numberOfTransitions(); i++)
 
427
                set_col_name(lp, i+1, const_cast<char*>(net.transitionNames()[i].c_str()));
 
428
 
 
429
        // Start adding rows
 
430
        set_add_rowmode(lp, TRUE);
 
431
 
 
432
        REAL row[net.numberOfTransitions() + 1];
 
433
        for(size_t p = 0; p < nPlaces; p++){
 
434
                // Set row zero
 
435
                memset(row, 0, sizeof(REAL) * net.numberOfTransitions() + 1);
 
436
                for(size_t t = 0; t < net.numberOfTransitions(); t++){
 
437
                        int d = net.outArc(t, p) - net.inArc(p, t);
 
438
                        row[1+t] = d;
 
439
                }
 
440
 
 
441
                if(pcs[p].min == pcs[p].max &&
 
442
                   pcs[p].max != CONSTRAINT_INFTY){
 
443
                        double target = pcs[p].min - m0[p];
 
444
                        add_constraint(lp, row, EQ,  target);
 
445
                }else{
 
446
                        // There's always a min, even zero is interesting
 
447
                        double target = pcs[p].min - m0[p];
 
448
                        add_constraint(lp, row, GE,  target);
 
449
                        if(pcs[p].max != CONSTRAINT_INFTY){
 
450
                                double target = pcs[p].max - m0[p];
 
451
                                add_constraint(lp, row, LE,  target);
 
452
                        }
 
453
                }
 
454
        }
 
455
 
 
456
        // Finished adding rows
 
457
        set_add_rowmode(lp, FALSE);
 
458
 
 
459
        // Create objective
 
460
        memset(row, 0, sizeof(REAL) * net.numberOfTransitions() + 1);
 
461
        for(size_t t = 0; t < net.numberOfTransitions(); t++)
 
462
                row[1+t] = 1;   // The sum the components in the firing vector
 
463
 
 
464
        // Set objective
 
465
        set_obj_fn(lp, row);
 
466
 
 
467
        // Minimize the objective
 
468
        set_minim(lp);
 
469
 
 
470
        // Set variables as integer variables
 
471
        for(size_t i = 0; i < net.numberOfTransitions(); i++)
 
472
                set_int(lp, 1+i, TRUE);
 
473
 
 
474
        // Attempt to solve the problem
 
475
        int result = solve(lp);
 
476
 
 
477
        // Limit on traps to test
 
478
        size_t traplimit = nPlaces * OVER_APPROX_TRAP_FACTOR;
 
479
        // Try to add a minimal trap constraint
 
480
        while((result == OPTIMAL) && traplimit-- < 0){
 
481
                memset(row, 0, sizeof(REAL) * net.numberOfTransitions() + 1);
 
482
                // Get the firing vector
 
483
                get_variables(lp, row);
 
484
                // Compute the resulting marking
 
485
                MarkVal rMark[net.numberOfPlaces()];
 
486
                for(size_t p = 0; p < nPlaces; p++){
 
487
                        rMark[p] = m0[p];
 
488
                        for(size_t t = 0; t < net.numberOfTransitions(); t++)
 
489
                                rMark[p] += (net.outArc(t, p) - net.inArc(p, t)) * (int)row[t];
 
490
                }
 
491
 
 
492
                // Find an M-trap
 
493
                BitField trap(minimalTrap(net, m0, rMark));
 
494
 
 
495
                //Break if there's no trap
 
496
                if(trap.none()) break;
 
497
 
 
498
                // Compute the new equation
 
499
                for(size_t t = 0; t < net.numberOfTransitions(); t++){
 
500
                        row[1+t] = 0;
 
501
                        for(size_t p = 0; p < nPlaces; p++)
 
502
                                if(trap.test(p))
 
503
                                        row[1+t] += net.outArc(t, p) - net.inArc(p, t);
 
504
                }
 
505
 
 
506
                // Add a new row with target as greater than equal to 1
 
507
                set_add_rowmode(lp, TRUE);
 
508
                add_constraint(lp, row, GE,  1);
 
509
                set_add_rowmode(lp, FALSE);
 
510
 
 
511
                // Attempt to solve the again
 
512
                result = solve(lp);
 
513
        }
 
514
 
 
515
        int retval = 0;
 
516
 
 
517
        if(result != INFEASIBLE){
 
518
                get_variables(lp, row);
 
519
                for(size_t t = 0; t < net.numberOfTransitions(); t++)
 
520
                        retval += (int)row[t];
 
521
        }
 
522
 
 
523
        // Delete the linear problem
 
524
        delete_lp(lp);
 
525
        lp = NULL;
 
526
 
 
527
        // Return true, if it was infeasible
 
528
        return retval;
 
529
}
 
530
 
 
531
} // Structures
 
532
} // PetriEngines